Dark matter in the solar system I: The distribution function of WIMPs at the Earth 

from solar capture. 
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The next generation of darlc matter (DM) direct detection experiments and neutrino telescopes 
will probe large swaths of dark matter parameter space. In order to interpret the signals in these 
experiments, it is necessary to have good models of both the halo DM streaming through the solar 
system and the population of DM bound to the solar system. In this paper, the first in a series of 
three on DM in the solar system, we present simulations of orbits of DM bound to the solar system 
by solar capture in a toy solar system consisting of only the Sun and Jupiter, assuming that DM 
consists of a single species of weakly interacting massive particle (WIMP). We describe how the size 
of the bound WIMP population depends on the WIMP mass m^, spin- independent cross section 



Up , and spin-dependent cross section . Using a standard description of the Galactic DM halo, 
we find that the maximum enhancement to the direct detection event rate, consistent with current 
experimental constraints on the WIMP-nucleon cross section, is < 1% relative to the event rate 
from halo WIMPs, while the event rate from neutrinos from WIMP annihilation in the center of the 
Earth is unlikely to meet the threshold of next-generation, km'^-sized (IceCube, KMSNeT) neutrino 
telescopes. 

PACS numbers: 95.35.+d,96.25.De,95.85.Ry,96.60.Vg 



I. INTRODUCTION 



Dark Matter and Detection 



There is overwhelming evidence that non-baryonic DM 
must exist in large quantities in the universe, yet its na- 
ture is unknown. A popular candidate for DM is one or 
more species of WIMP. Particles of this type occur nat- 
urally in many theories of physics beyond the Standard 
Model (SM); examples include the neutralino x iii super- 
symmetry [1], the lightest Kaluza-Klein photon B^'^^ in 
universal extra-dimension (UED) theories [2-4], or the 
heavy photon Ah in Little Higgs models [5-7]. These 
particular candidates are all stable, self-annihilating, be- 
have as cold dark matter, and are thermally produced 
in the early universe in roughly the amount needed to 
explain the dark matter [8]. 

We may expect rapid progress in constraining the na- 
ture of DM due to the maturity of a number of tech- 
nologies targeting different but complementary WIMP 
signals. The next generation of particle colliders, in par- 
ticular the Large Hadron Collider, may see signatures of 
physics beyond the SM. A new generation of satellites is 
searching for photons (e.g., the Fermi Gamma-Ray Space 
Telescope [9, 10]) and other particles (e.g., ATIC [11], 
PAMELA [12]) resulting from WIMP annihilations in the 
Milky Way's DM halo. 

There are also experiments to probe the local WIMP 
population. Since the flux of particles from WIMP an- 
nihilation scales as the square of the WIMP density, any 



region in the solar system that has an unusually high 
density of WIMPs is a good target. WIMPs generically 
interact with baryons, which means that WIMPs pass- 
ing through the solar system may be trapped and settle 
into dense cores in the potential wells of the Sun or the 
planets. The previous generation of neutrino telescopes 
(e.g., BAKSAN [13], MACRO [14], Super-Kamiokande 
[15], AMANDA [16, 17]) places the strongest constraints 
on the spin-dependent WIMP-proton cross section 
(< 10"'^^ cm^ at ^ 100 GeV) based on flux limits of 
neutrinos from WIMP annihilation in the Sun and Earth. 
Even if the next generation of neutrino telescopes with 
detector volumes approaching 1 km'^ (e.g., Antares [18], 
IceCube [19], the proposed KMSNeT [20]) do not iden- 
tify a WIMP annihilation signature, they are projected 
to improve constraints on a^^ by almost two orders of 
magnitude. 

The best limits on the spin-dependent WIMP-neutron 
and spin-independent WIMP-nucleon a^^ cross sec- 
tions come from direct detection experiments. The sig- 
nature of WIMPs in these experiments is a small (~ 
10 keV — 100 kev) nuclear recoil. The next generation 
of direct detection experiments is slated to have target 
masses approaching 1000 kg (e.g., DEAP/CLEAN [21], 
LUX [22], SuperCDMS [23-25], WARP [26], XENONIT 
[27], XMASS [28]) and to be sensitive to cross sections 
down to ^ cm^, cr^^ - lO""*" cm^, and 
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1^ cm"^ 



[29, 30]. 



B. WIMPs in the Solar System 



•Electronic address: apeteraastro.caltech.edu 



For a given WIMP model, event rates in direct detec- 
tion experiments and neutrino telescopes are determined 
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by the phase space distribution function (DF) of WIMPs 
in the solar system. The fiducial assumption is that the 
direct detection event rate and WIMP capture rate in the 
Earth arc dominated by DM particles from the Galactic 
halo, passing through the solar system on unbound orbits 
[1, 31]. There are potentially observable consequences if 
even a tiny fraction of WIMPs may become captured to 
the solar system, since bound WIMPs have lower speeds 
than halo WIMPs. The push for many direct detection 
experiments is toward ever-lower nuclear recoil energy 
thresholds (~ 0.1 keV — 1 keV), both in order to gain 
sensitivity to low mass WIMPs and because the event 
rate is much higher there than at lower energies [32, 33]. 
At such low energies, low speed WIMPs contribute dis- 
proportionately to the event rate for kinematic reasons. 

Low speed WIMPs have an even greater impact on the 
event rate of neutrinos from annihilation in the Earth. 
The shallowness of the Earth's potential well means that 
only low speed WIMPs may be captured in the Earth. 
In particular, if the WIMP mass is above 400 GeV, only 
WIMPs bound to the solar system may be trapped in the 
Earth. 

Two processes have been identified by which WIMPs 
may become captured to the solar system at rates large 
enough to be important for terrestrial dark matter exper- 
iments. Gravitational Capture: Gould [34, 35] pointed 
out that WIMPs may be captured from the halo by grav- 
itationally scattering on the planets. By treating WIMP 
orbits in the solar system as a diffusion problem, Gould 
[35] and Lundberg and Edsjo [36] estimated that bound 
WIMPs dominate the annihilation rate of WIMPs in the 
Earth for WIMP masses > 100 GcV. Solar Capture: 
WIMPs captured in the Sun will reach thermal equilib- 
rium with solar nuclei on timescales t t~^P, where 
T is the optical depth of the Sun for WIMPs and P is 
the orbital period of a bound WIMP. However, Damour 
and Krauss [37] identified a population of solar-captured 
WIMPs that could survive for much longer periods of 
time due to a type of secular resonance that pulls their 
perihelia outside the Sun. Using secular perturbation 
theory, they found that this population could produce a 
low-recoil direct detection rate comparable to that of halo 
WIMPs for <7^^ - 10-'*2 _ 10-40 cm2, and could yield an 
annihilation rate in the Earth a factor of ~ 100 greater 
than the rate expected from unbound halo WIMPs for 
WIMP masses - 100 - 150 GeV [38]. 

While these results are intriguing, the semi-analytic 
treatments used in these papers cannot fully capture the 
rich range of behavior in small-N systems such as the so- 
lar system. It is important to check these results with nu- 
merical experiments. Moreover, the annihilation rate of 
WIMPs in the Sun depends critically on whether WIMPs 
captured in the Sun thermalize rapidly with solar nuclei. 
If the planets can pull the WIMPs out of the Sun for 
extended periods of time, or even eject the particles from 
the system, the annihilation rate will be depressed with 
respect to current estimates. 

In a set of three papers [39, 40], we present simulations 



of WIMP orbits in the solar system, including both the 
gravitational efi^ects of the dominant planet, Jupiter, and 
an accurate Monte Carlo description of WIMP-nuucleon 
elastic scattering in the solar interior, as well as a discus- 
sion of the likely contribution of bound WIMPs to direct 
detection experiments and neutrino telescopes. In this 
paper, Paper I, we focus on WIMP capture in the Sun. 
In order to put our results in context, we first summarize 
the treatment of Damour and Krauss [37] , and describe 
the mechanism they found that extended the lifetimes of 
solar captured WIMPs in the solar system and built up 
the DF of WIMPs at the Earth: the Kozai mechanism. 



C. Damour and Krauss (1999) and the Kozai 
Mechanism 

In the absence of gravitational torques from the plan- 
ets, WIMPs captured onto Earth-crossing orbits by elas- 
tic scattering in the Sun will have a small number density 
at the Earth relative to the halo number density for two 
reasons, (i) Unless the WIMP is massive (m > 1 TeV), 
the characteristic energy a WIMP loses to a solar nu- 
cleus is large enough such that most captured WIMPs 
have aphelia that lie inside the Earth's orbit, (ii) The 
characteristic time to the next scatter, which will almost 
certainly remove the WIMP from an Earth-crossing orbit 
unless m > 1 TcV, is of order t oc P/t. For a WIMP with 
semi-major axis a = 1 AU, P^ — 1 year. If, for example, 
a^i = 10-41 cm2 (or = lO'^s cm^), r ~ lO'^, so 
the WIMP lifetime in the solar system is only of order a 
thousand years, short compared to the age of the solar 
system. 

Damour and Krauss [37] recognized that the lifetimes 
of bound WIMPs in Earth-crossing orbits could be; ex- 
tended by orders of magnitude if gravitational torques 
from the planets decreased the WIMP eccentricity (in- 
creased the perihelion distance) enough that the WIMP 
orbit no longer penetrated the Sun. For WIMPs in plan- 
etary systems such as our own, such behavior is possible 
if the rate of perihelion precession lu is small, since then 
the torques from the planets act in a constant direction 
over many WIMP orbits. This process was first exam- 
ined by Kozai [41] in the context of asteroid orbits and 
is sometimes called the Kozai resonance. The signature 
of this resonance is large fluctuations in both the inclina- 
tion and eccentricity while the semi-major axis is fixed. 
The Kozai resonance can lead to both libration and cir- 
culation in the argument of perihelion lu, and we use the 
term "Kozai cycles" to describe these oscillations. Kozai 
cycles have been studied in the context of comets [42], as- 
teroids [43, 44], triple star systems [45], and exoplanets 
[46]. 

Damour & Krauss found approximate analytic Kozai 
solutions for a solar system containing the inner planets 
and Jupiter on circular, coplanar orbits. The require- 
ment that u! is small means that Kozai cycles are only 
significant for WIMPs with perihelia that are not too 
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far inside the solar radius, so the solar potential is not 
far from that of a point mass. Damour & Krauss use 
an analytic approximation to the solar potential in the 
outer r > 0.55Rq of the Sun, where Rq is the radius 
of the Sun. They expanded the potentials of the plan- 
ets to quadrupole order in the small parameter ap/a, 
where ap is the semi- major axis of a planet, and ne- 
glected short-period terms and mean-motion resonances. 
The solutions have an additional feature — if the orbital 
plane of the planets in the solar system is the x — y 
plane, the 2:- component of t he specific angular momen- 
tum, Jz = y^GMQa{l — e^) cos / (7 is the inclination) is 
a conserved quantity. 

To estimate the size of the solar-captured WIMP pop- 
ulation at the Earth, Damour & Krauss made the follow- 
ing additional assumptions, (i) Jupiter-crossing WIMPs 
(with aphelia greater than Jupiter's semi-major axis, 
fa > Qi). ~ 5.2 AU) were ignored, since their lifetimes 
in the solar system were assumed to be short. Similarly, 
all WIMPs with < aq^ that were not on Kozai cycles 
were also ignored, (ii) They assumed that all WIMPs 
on Kozai cycles would survive for the lifetime of the so- 
lar system without rescattering in the Sun, regardless of 
the optical depth in the Sun for WIMPs. Since the typ- 
ical lifetime of Earth-crossing WIMPs is ~ 10^ yr for 
^ 10~^^ cm^, the extension of the lifetimes of even a 
few Earth-crossing particles to the age of the solar system 
results in a significant boost to the DF of bound WIMPs 
at the Earth. 



D. This Work 

To investigate the validity of these assumptions and 
to provide a more accurate assessment of the contribu- 
tion of bound WIMPs to direct detection experiments 
and neutrino telescopes, we perform a set of numerical 
simulations of WIMP orbits in the solar system. In this 
paper. Paper I, we present a suite of simulations of solar- 
captured WIMP orbits in a toy solar system consisting 
only of the Sun and Jupiter. Jupiter is the only planet 
included in the simulations for two reasons, (i) As the 
largest planet in the solar system, it dominates the dy- 
namics of minor bodies in the system. We address the 
issue of other planets in Section VII of this paper as 
well as in later papers, (ii) Since some of our numeri- 
cal methods (described in Section II) are new, and since 
it is important to have a physical understanding of the 
principal mechanisms that determine the key features of 
the bound WIMP population, it is useful to simulate a 
simple system first. In particular, particle orbits in our 
toy solar system enjoy a constant of motion (Eq. 20), 
which provides a check on the numerical accuracy of the 
integrations. 

We describe the simulations in Section III, and the 
DFs derived from the simulations in Section IV. Also in 
Section IV, we show how the DFs depend on the WIMP 
mass and cross section m^, cTp^, and Up^. Predictions 



for the event rates in direct detection experiments and 
neutrino telescopes are made in Sections V and VI. We 
discuss our results in the context of the previous work on 
solar captured- WIMPs by Damour and Krauss [37] and 
Bergstrom et al. [38] , the presence of other planets, and 
the assumptions concerning the halo DF in Section VII, 
and summarize the main results of this work in Section 
VIII. 

We defer the topic of annihilation of WIMPs in the Sun 
to Paper II [39], and the simulations of gravitationally 
captured WIMPs to Paper III [40]. 



II. ORBIT INTEGRATION 

The problem of determining the long-term trajectories 
of bound dark matter particles imposes a set of difficult 
challenges to the integration algorithm. The algorithm 
must (i) be stable and accurate over 4.5 Gyr; (ii) accu- 
rately follow highly eccentric (e > 0.995) orbits with no 
numerical dissipation; (iii) accurately integrate trajecto- 
ries that are influenced by perturbing forces that may be 
comparable to the force from the Sun for short intervals 
(including close encounters with and passages through 
planets); and (iv) be fast, in order to obtain an adequate 
statistical sample of orbits. 

Much progress has been made in the past fifteen years 
to address the first and last criteria. This progress has 
largely been motivated by interest in the long-term sta- 
bility of planetary systems. The most significant develop- 
ment has been the advent of geometric integrators (sym- 
plectic and/or time-reversible integrators), which have 
the desirable property that errors in conserved quanti- 
ties (such as the Hamiltonian) are oscillatory rather than 
growing. However, the most commonly used algorithms 
[47-49] are not immediately applicable to the present 
problem, for two main reasons. First, one would like 
to use an adaptive time step to quickly but accurately 
integrate a highly eccentric orbit (using very small time 
steps near perihelion and larger ones otherwise), or to re- 
solve close encounters with the planets. It is difficult to 
introduce an adaptive time step in a symplectic or time- 
reversible way since varying the time step by criteria that 
depend on phase space position destroys symplecticity. 
Secondly, since for practical purposes the integrations of 
planetary or comet orbits end when two bodies collide, 
there has been little attention to integrating systems for 
which the potential can deviate significantly from the Ke- 
plerian point-mass potential, as it does in the solar inte- 
rior. 

In the following sections, we describe an algorithm to 

efficiently carry out the long-term integration of dark 
matter particles in the solar system. In Section II A, 
we outline an adaptive time step symplectic integrator 
(simultaneously formulated by Preto and Tremaine [50] 
and Mikkola and Tanikawa [51]) that is used for most 
of the orbital integrations. We explain the error proper- 
ties of the integrator in Section II B. In Section II C, we 
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discuss procedures to handle special cases, such as close 
planetary encounters. We discuss the merits of various 
coordinate systems in Section II D. 



A. The Adaptive Time Step Integrator 

We closely follow the arguments of Mikkola and 
Tanikawa [51] and Preto and Tremaine [50] in the de- 
scription of the adaptive time step symplectic integrator. 

A separable Hamiltonian -ff(q, p, t) = T{p) + U (q, t) 
(T is the kinetic energy and U is the potential energy), 
a function of the canonical position q and momentum p, 
can be implemented as a symplectic integrator with fixed 
time step At. The key to finding a symplectic integrator 
with a variable time step is to promote the time t to 
a canonical variable, and make it a function of a new 
"time" coordinate s, 



At = ,g(q,p,t)ds. 



(1) 



find a separable Hamiltonian F in the extended phase 
space that describes the motion, and take fixed time steps 
As when integrating the new equations of motion. The 
new canonical coordinates are qo = t and po = —H, so 
the new set of canonical variables Q = {qo,o^ and P = 
(pojP)- Preto and Tremaine and Mikkola and Tanikawa 
find such an extended phase space Hamiltonian, 



r(Q,P)=5(Q,P) [//(q,P,i)-po], 
which can be made separable with the choice 

- /(np)+Po)-/(-^(Q)) 
^^""^'"^^^ T(p) + [/(Q)+po ' 

so that the extended Hamiltonian is 

r(Q,p) = /(r(p)+po)-/(-c/(Q)). 

The equations of motion for this Hamiltonian are: 
/'(T(p)+po) 



dgo 
ds 
dq 

ds 
dpo 
ds 
dp 
ds 



/'mP)+Po)^ 



-/'(-c/(Q) 



dqo 
9q 



(2) 
(3) 

(4) 

(5) 
(6) 
(7) 
(8) 



To determine a useful choice for f{x), Preto and 
Tremaine expand Eq. (3) in a Taylor scries about the 
small parameter T + po + U (= if the Hamiltonian is 
exactly conserved) to show that 



Outside the Sun, the gravitational potential of the solar 
system is 



C/(q,t) 







jq-q© 



(10) 



where the first term in the potential denotes the Kep- 
lerian potential of the Sun and is the potential from 

planet i, and the potential from the Sim dominates most 
of the time. Preto and Tremaine show that for a choice 
of 



9(Q,P) 



[q-q©! 

GMq 



(11) 

(12) 



the two-body problem can be solved exactly, with only a 
time (phase) error St/P oc A''"^, where P is the orbital 
period and N is the number of steps per orbit. This is a 
good feature because phase errors are far less important 
for our purposes than, for example, systematic energy 
drifts or numerical precession. Note that the time step 
is proportional to the particle's separation from the Sun, 
so that small time steps are taken near the perihelion of 
the orbit and large steps near the aphelion. We use Eq. 
(11) as our choice for g{Q,P), for which the functional 
form of f{x) is 



f{x)=GMQlog{x). 



(13) 



The adaptive time step equations of motion are imple- 
mented via a second-order leapfrog integrator (also called 
a Verlet integrator) with As ~ At/g = /i, where h is de- 
termined by the number of steps desired per orbit. Since 
the goal is to understand the behavior of a large ensem- 
ble of orbits, we are more interested in maintaining small 
energy errors over long times rather than precisely in- 
tegrating orbits over short times, and so a second-order 
symplectic integrator is sufficient. For our choice of f{x), 
and given T = w^/2 and U — U{r,t), the change over a 
single fictitious time step h can be written as the map- 
ping 



ri/2 

tl/2 



Vl = 



POA 



ri = 



1 GMqVo 

2 2^0+^0,0 
1, GMq 

2 2^0 + Po,o 

vo + h— — ^ ^ '— 

GMq dU{Vi,2,h/2) 

Po,o + h— — -7 

U{ri/2,ti/2) at 

''1/2+ 



tl/2 



2^1 +P0,1 
1, GMq 

2 H+Po,i' 



(14) 
(15) 
(16) 
(17) 
(18) 
(19) 



5(Q,P)«/'(-c/(Q))- 



(9) 



where the subscript i = 0, 1/2, 1 labels what fraction of 
the total time step h has been taken. 
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B. Errors Along the Path 



We explore the behavior of the energy errors in the 
adaptive time step integrator as a function of energy, 
eccentricity, distance from the Sun, and number of steps 
per orbit. This study allows us, in conjunction with the 
results of Section II C 2, to determine which the fictitious 
time step h to use to meet accuracy requirements. The 
choice of h for the simulations is described in Section 
HID. For the current study, we use short integrations 
in order to focus on the errors of the adaptive time step 
integrator alone. We will discuss the long-term behavior 
of the whole integration scheme after we discuss the other 
pieces of our algorithm. 

Since our toy solar system (Sun + Jupiter + WIMP) 
is a restricted three-body problem, there is one constant 
of motion, the Jacobi constant 



Cj = -2(S-n^J,), 



(20) 



where E is the particle energy in an incrtial frame, nq^ is 
the mean motion of Jupiter, and Jz is the ^-component 
of the particle's angular momentum, assuming that the 
motions of the Sun and Jupiter are confined to the x — y 
plane. Therefore, we can parameterize errors in terms of 
the Jacobi constant. There are no analogous conserved 
quantities for particles orbiting in planetary systems with 
more than one planet or if the planetary orbit is eccentric. 
In those systems, one can quantify errors for integrators 
of the type described in Section II A in terms of the rela- 
tive energy error AE/E = {E+po)/E, where E is deter- 
mined by the instantaneous position and velocity of the 
particle and po is the 0— component of the momentum 
in the extended phase space. If the equations of motion 
(5)-(8) were integrated with no error, then po — —E and 
AE/E = 0. 

In this section, we treat the Sun as a point mass, and 
consider trajectories with aphelia well inside Jupiter's or- 
bit. We consider two different initial semi-major axes, 
a = a7j_/3 and a = aq^/G respectively, where otj, is the 
semi-major axis of Jupiter. To determine the size of the 
errors in Cj as a function of eccentricity, we integrate 
orbits with initial eccentricity e = 0.9, 0.99, 0.999 and 
0.9999. We perform integrations for each combination of 
a and e for 10 different initial, random orientations and 
an ensemble of step sizes. We run each integration for 
a total of 2 X 10** Kepler periods. The integrations are 
started at perihelion (to mimic the initial conditions of 
scattering in the Sun) with a very small h — IO^^Rq^ 
year. We use such a small time step because the mag- 
nitude of the errors in the integrator are largest if the 
integration is started at pericenter, and smallest when 
started at apocenter. Once the particle reaches its first 
aphelion, h is adjusted so that it will provide the desired 
number of steps per orbit. The fictitious time step is re- 
lated to the number of steps per orbit by the step in the 
eccentric anomaly Ait and semi-major axis a by 

1 — cos Am 
(GMq /a) 1/2 sin Au' ^ ^ 




FIG. 1: Jacobi constant errors as a function of distance from 
the primary for a trajectory with a = 1.73 AU, followed for 
2 X 10* Kepler periods. This trajectory was integrated with 
500 steps/orbit. Errors are calculated at perihelion and aphe- 
lion. Points to the left of the vertical line lie within the vol- 
ume of the Sun; however, we used a point-mass Sun for this 
integration. 



for the symplectic mapping of Eqs. (14)-(19) in the case 
of the Kepler two-body problem. The number of steps 
per orbit is given by 



N = 



2n 
Au' 



(22) 



We show the dependence of the error on the distance 
from the Sun in Fig. 1. In this figure, we plot the perihe- 
lion and aphelion Jacobi constant errors for a trajectory 
with initial a — aq^/3 and e = 0.999, integrated with 
500 steps/orbit, representative of all the simulations. We 
plot only errors at perihelion and aphelion for clarity; a 
plot showing errors at each time step would be similar 
but with more scatter. The interior of the Sun is in the 
shaded region (though the integrations were done for a 
point-mass Sun). From Figure 1, it appears that 



\ACj/Cj\^r' 



(23) 



This is a generic feature of the integrator, and implies 
that the maximum Jacobi constant or energy error occurs 
at perihelion. The errors are oscillatory, i.e., there is no 
secular growth in the error envelope with time. 

In Fig. 2, we show the maximum Jacobi constant error 
as a function of initial semi-major axis and eccentricity 
Ci. To find the maximum error, we calculate the error in 
Jacobi constant every time e is in the range €^±0.1(1— e^). 
The restriction on e isolates the effect of eccentricity on 
\ACj /Cj\, since Fig. 1 demonstrates that the maximum 
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FIG. 2: Errors in the Jacobi constant as a function of eccen- 
tricity and semi-major axis. Each point shows the maximum 
error for 10 trajectories initialized with the same eccentricity 
but with random initial orientation, and followed for 2 x 10'* 
Kepler orbits. Open points denote those trajectories for which 
the semi-major axis a = a'^_/3 = 1.73 AU; closed points refer 
to trajectories with a = aq^./6 = 0.87 AU. Circles mark tra- 
jectories with initial eccentricity Ci = 0.9999, squares denote 
those with a = 0.999, diamonds indicate those with a = 0.99, 
and triangles those with a = 0.9. 



error in a simulation depends on the largest eccentricity 
in the orbit. Wc then plot the maximum error found 
among all simulations for the same initial and e^. For 
each type of simulation, the maximum error occurs at 
perihelion. Fig. 2 indicates that the maximum Jacobi 
constant error is a decreasing function of the number of 
steps per orbit, and an increasing function of semi-major 
axis and eccentricity. Furthermore, the maximum error 
for e S ej±0.1(l — Ci) within each simulation is a function 
of the initial conditions. In the simulations with fixed 
eccentricity and a = aq^_/6, the spread in these central 
values is less than a factor of two, while the spread is 
about a factor of ten in the a = ai^/S simulations. This 
is described more in [52]. 

To set the fictitious time step h for the simulations de- 
tailed in Section HID, it is preferable to consider errors 
at a fixed, small distance from the Sun rather than ex- 
clusively at perihelion. This is because we use a mapping 
technique to follow perihelion passages where r-p < 2Rq. 
Therefore, we want to impose a maximum Jacobi con- 
stant (or energy) error for the simulations at r = 2Rq. 
However, we also want to optimize h such that passages 
near planets can be integrated accurately with the least 
overall CPU time. A full discussion which values h are 
used for the main set of simulations in this work will be 



deferred to Section HID, after we discuss close encoun- 
ters with Jupiter in Section II C 2. 



C. Special Cases 

While we would like to use this adaptive time step 
integrator as much as possible, keeping the fictitious step 
h fixed, there are two situations which must be handled 
separately. 



1. The Sun 

The interior of the Sun has a potential that deviates 
strongly from Keplerian. The integrator described in Sec- 
tion II A works badly inside the Sun because it is designed 
for nearly Keplerian potentials. Thus, we replace the in- 
tegration through the Sun by a map. We exploit the fact 
that tidal forces from the planets are small near the Sun. 
Since the two-body problem can be solved exactly, we can 
define a region about the Sun (called a "bubble," with 
a typical radius of 0.1 AU) for which we use the exact 
solutions to the two-body problem. In reality, we create 
a map for the bubble but only use it if the orbital perihe- 
lion lies within r = 2Rq. The bubble wall is larger than 
2Rq so that a particle does not accidentally step into the 
Sun when stepping into the bubble. In the WIMP orbital 
plane, we map the incoming position and velocity to the 
outgoing position and velocity using look-up tables for 



At{a, e) = 



rn 
J Tp{a,e 



dr 



e) J2[±l-~^Q{T)]±a{l-e^)/T^ 



(24) 



and 



A0(a, e) = 2^±a{e^ - 1) 

Jrp{a,e) / 



dr 



2f±. 



$0(r)]±a(l-e2)/r2 



(25) 



which are the time At and phase Acp through which the 
particle passes in the bubble region. By convention, a is 
always positive, such that E = GMQ/2a for hyperbolic 
orbits and E = —GMQ/2a for eccentric orbits. The +/— 
signs in Eqs. (24) and (25) correspond to hyperbolic 
(e > 1) and elliptical orbits (e < 1) respectively. We have 
normalized the solar potential $0 = ^q/GMq. Note 
that Th is the bubble radius and Vp is the true perihelion, 
defined by 



dr 
dt 



= (26) 



/2(±--$0(r^))±a(l-e2)/r2.(27) 
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We parameterize the look-up tables in terms of the semi- 
major axis and Keplerian perihelion vk = — 

There is one subtlety in matching the map through 
the bubble to the integrator outside of bubble. In the 
Keplerian two-body problem, one solves the equations of 
motion dp/dt and dq/dt instead of dp/ds and dq/ds. If 
one divides dq/ d.s and dp / ds by the differential equation 
for the time coordinate, the time-transformed equations 
of motion are 



dq 

dt 



dp 

dt 



dq/ds 
dt/ds 

f'{T + po)dT/dp 
f'iT + Po) 



dp/ds 

dt/ds 

f{-U)dU/dci 

f'{T + Po) 
f{-U) dU 
/'(T + po)"9q' 



(28) 

(29) 

(30) 

(31) 

(32) 

(33) 



The second of these differs from the equations of motion 
of the original Hamiltonian by a multiplicative factor 



^l = .n~u)/.f'{T + po), 



(34) 



in other words, the particle follows a Kepler orbit about 
a Sun of mass iiMq. Therefore, we calculate the orbital 
elements using 



a = 



Po 



e = ^/lny{jIGM^, 



(35) 

(36) 



where the upper (lower) sign should be used for hyper- 
bolic (elliptical) orbits. We use a look-up table for At 
and A(j) with the modification that At, as calculated for 
a and e with n = 1, must be multiplied by a factor of 
^-i/2_ rjnj^^ change in phase is unaffected by the choice 
of central mass since it is a purely geometric quantity. 

Similar lookup tables are also used to determine the 
perihelion and the optical depth as a function of semi- 
major axis and eccentricity. We discuss additional scat- 
tering in the Sun in Appendix B. 

We demonstrate the robustness of the map in the up- 
per left panel of Fig. 3, where we show errors in the 
Jacobi constant over a 500 Myr time span for an orbit 
with a ~ 1.54 AU. The orbit enters the Sun ~ 10'' times 
in this time span. We sample the orbit at the first aphe- 
lion after a 10^ yr interval from the previous sample, 
and there are approximately 100 steps/orbit. This fig- 
ure shows that there are only oscillatory errors through- 
out this long-term integration, and these fractional errors 
never exceed 10~^ at aphelion. Long-term integrations of 
the two-body problem using the map demonstrate energy 
errors only at the level of machine precision. 



2. The Planets 

While the adaptive time step integrator works well in a 
near-Keplerian potential, one must treat close encounters 
with planets more carefully. If the time step is too large 
near a planet, the particle fails to resolve the force from 
the planet. This can cause growing errors in the parti- 
cle's trajectory. Since we use an f{x) that is optimized to 
the potential of the Sun, the only way to achieve a small 
time step near each planet is to cither make the fictitious 
time step h small or to switch to a different integration 
method near each planet while using the method of Sec- 
tion II A with a reasonably large h for the rest of the 
orbit. The advantage of the former approach is that it 
does not break the symplectic nature of the integrator. 
However, it is also prohibitively computationally expen- 
sive. Therefore, we use the latter approach. 

Wc define a spherical region ("bubble") about each 
planet for which we allow a different integration scheme, 
while continuing to use the adaptive time step symplec- 
tic integrator (Section II A) outside the spheres. The 
transition between the integration schemes is not sym- 
plectic, but reduce errors in the integration by enforcing 
an accuracy requirement on \AE/E\ = \{pq + E)/E\ = 
\{—H + E)/E\ in the bubble of each planet. 

In the bubble of each planet, we continue to use the 
adaptive time step integrator, but the value of h' (the 
prime denotes the fact that this fictitious time step is 
only used within a planet bubble) used in the bubble is 
selected to keep the quantity \ AE / E\ as small as possible 
while also keeping the total integration time short. To 
find the optimal value of h' , we use the following algo- 
rithm. When a particle first enters a bubble, we record 
the particle's energy error at the first step, \AEi/Ei\. 
Then, we integrate the particle's trajectory through the 
bubble using the default value oi h. As the particle is 
about to exit the bubble, we calculate the energy error 
\AE f / Ef\. If the energy error meets the accuracy crite- 
rion, or if it is less than \AEi/Ei\, then the integration 
is allowed to continue normally. If, however, |A£'//£'/| 
does not satisfy the accuracy criterion, we restart the 
integration in the bubble from the point at which the 
particle first entered with a smaller fictitious time step 
h' . This process iterates until either the energy accu- 
racy condition is satisfied or the energy error plateaus 
in value. If the energy error plateaus in value, whichever 
trajectory (corresponding to a particular choice of h') has 
the minimum \AEf /Ef \ is chosen. 

The choice of the bubble size 1% is related to the choice 
of fiducial value of h and to the mass of the planet. A 
larger value of h means that the bubble needs to be larger 
to ensure that the planet's gravitational potential is prop- 
erly resolved. Planets with larger masses will require 
larger bubbles than smaller planets. We choose to keep 
the bubble size fixed for all orbits. In general, we tune 
h so that the typical energy errors for all energies are 
similar near each planet, and to keep the error in the Ja- 
cobi constant small \ACj/Cj\ < lO"" at r = 2Rq. The 
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FIG. 3: Error in the Jacobi constant as a function of time for several particles. The Jacobi constant is recorded at aphelion 
at 10^ yr intervals. Top left: A particle with a = 1.54 AU. This particle repeatedly goes through the Sun (about 10^ times), 
but never goes through the bubble around Jupiter. It is integrated with ^ = 6 x 10~*-Rq^ yr, which corresponds to ~ 100 
steps/orbit. Top right: A particle that gets stuck near a Sun-skimming 2:1 resonance with Jupiter. This particle repeatedly 
goes through the Jupiter bubble. It is integrated with h — 2 x 10~^-Rq^ yr, or ~ 350 steps/orbit. Bottom left: A particle gets 
stuck near a 3:2 resonance with Jupiter. This orbit was integrated with h — 1.5 x 10~^f?g^ yr, or ~ 650 steps/orbit. Bottom 
right: This particle repeatedly crosses rc, the transition radius between barycentric and heliocentric coordinates (dashed line 
marks rc/2, the crossing semi-major axis for an orbit with e ~ 1) and has its last aphelion before ejection from the solar system 
at t = 1.6 X 10® years. It is integrated with h = 2 x yr, or 9 x 10^ steps/orbit. 



optimum sizes of the Jupiter bubble is ^ 1 — 3 AU if 
we require that \AEf/Ef \ - 10"'^ - 10"^ 

A complication arises when particles experience large 
changes in energy in their passage through the planetary 
bubble. In this case, the value of h that guaranteed a cer- 
tain precision in \AE/E\ in the pre-encounter orbit may 
be either too large (for adequate precision) or too small 
(it will slow down the orbital integration). Therefore, 
we change the value of h at the next aphelion. Again, 
this procedure breaks the symplectic nature of the inte- 
grator, but by changing h at aphelion, our experiments 
show that we minimize errors. In Section HID, we out- 
line how h is chosen for the initial orbits, and how h is 
changed if the particle experiences significant changes in 
energy from planetary encounters. 

We demonstrate the performance of the bubble for the 
case of the 3-body problem in Fig. 3. In this figure, the 
fractional error of the Jacobi constant is plotted against 
the time since the initial scatter in the Sun that pro- 
duced a bound orbit, and we show the first 500 Myr of 



the integrations. The Jacobi constant is measured at the 
aphelion of the orbit at lO'^ year intervals. The trajecto- 
ries of the particles in the upper right and bottom panels 
repeatedly pass through the bubble around Jupiter. For 
these integrations, 1% — 2.3 AU, and the energy criterion 
was \ AEf/Ef \ < 2 X 10~^. There are no secular changes 
of the Jacobi constant with time. Therefore, even though 
the planet bubble disrupts the symplecticity of the inte- 
grator, the integrator tracks the Hamiltonian well. 



D. Coordinate Choice 

For most of the integration, we use a heliocentric co- 
ordinate system for both the particles and the planets. 
There are two main reasons why we choose a heliocentric 
system. First, it is much simpler to use heliocentric coor- 
dinates for passages through the Sun. Secondly, consider 
the gravitational potential of the planets in the heliocen- 
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trie frame, 
$(r)p 



$d(r) + $i(r) 

GMp 

p 



\r-rp\ 



(37) 

-E^^, (38) 



where the indirect term («) arises from the fact that this 
coordinate system is not the center-of-mass coordinate 
system, and d denotes the direct term. For orbits that 
are well within a planet's orbit, the direct term can be 
expanded into spherical harmonics 



Mr) = E 

p 

= E 



GMp 



-rp\ 

GMp 
rp 



GMf 



GMf 
rp 



E 

1=2 



rp 



Pi 



-r • rp 



r • rp 



rrp 



(39) 
(40) 
(41) 



where the Pi are Legendre polynomials. The dipole term 
of the direct potential is canceled by the indirect poten- 
tial. Therefore, the primary contributor to the force on 
the particle by the planet comes from the I = 2 tidal term 
of the potential, whereas the I = 1 term is dominant in 
the barycentric frame. 

While there are distinct advantages to using the he- 
liocentric frame, the indirect term of the potential domi- 
nates the potential at large distances from the Sun. This 
poses a problem for the adaptive time step integrator, 
since the choice of g = —GMq/U = |r — r©! is only opti- 
mal if the Keplerian solar potential is dominant. There- 
fore, we choose to work in the barycentric frame at large 
distances. 

In practice, this means switching between heliocentric 
and barycentric coordinate systems for long-period or- 
bits. We choose the crossover radius such that 



max|$i,p(rc,^p = 0)| = e 



(42) 



where 9p is the angle between r and rp, rc is the 
crossover radius, the "max" signifies the planet for the 
planet P for which the indirect potential is strongest, 
and e is a factor < 1. In our solar system, the planet for 
which the indirect potential is strongest is Jupiter. The 
choice of e « 0.1 works well. The crossover radius is thus 



Mr. 



= e- 



or 



^eMQ/Mr^aq^. 



(43) 



(44) 



In changing coordinates, one breaks the symplectic flow 
of the integrator. Therefore, one must treat the Hamil- 
tonian, and therefore pq, carefully at the crossover. We 



choose to treat the transition the same way we treat the 
transition into the bubble about the Sun. Namely, we cal- 
culate (1 (Eq. 34), the factor by which the gravitational 
potential is modified in the integrator (see Eq. 33), in 
the initial coordinate frame i. Then we set 



Pol/ = -tJi.iE{x,t)\f, 



(45) 



where quantities calculated in the final frame are denoted 
by /, and E is the energy derived from the position and 
velocity coordinates of the particle. While this transition 
is not symplectic, in practice it conserves the Jacobi con- 
stant to adequate precision. This is demonstrated in the 
lower left panel in Figure 3, an orbit for which the initial 
semi-major axis is 50 AU. In this integration, e = 0.1, 
which translates to Tc = 53 AU. 



III. SIMULATIONS 
A. Dark Matter Model 

In order to perform the orbit simulations, it is neces- 
sary to specify some dark matter properties. The parti- 
cle mass and elastic scattering cross sections completely 
determine scattering properties in the Sun. and hence, 
these are the only WIMP-dependent parameters neces- 
sary to run the simulations and find the WIMP distribu- 
tion function at the Earth. The particle physics model 
and parameter space within each model do not need to 
be specified for the simulations, although we assume that 
the dark matter particle is a neutralino when we estimate 
event rates in neutrino telescopes in Section VI. Thus, 
we use to denote the WIMP mass. 

The relative strengths of the spin-dependent and spin- 
independent elastic scattering cross sections are impor- 
tant in the context of scattering in both the Sun and 
the Earth. For simplicity in interpreting the simulations, 
we would like to use either a spin-independent or spin- 
dependent cross section, but not a mixture of the two. 
We choose to focus on the spin-independent cross sec- 
tion for the simulations, but in Section IV D, we show 
how to extend our results to the case of non-zero spin- 
dependent interactions. In Section HID, we discuss the 
specific choices for the WIMP mass and a^^ used in the 
simulations. 

We adopt the Maxwellian distribution function (DF) 



v) 



(27ra2) 



^g-vV2<T'' 



(46) 



to describe the dark matter distribution function in the 

solar neighborhood in Galactoccntric coordinates and far 
outside the gravitational sphere of influence of the Sun. 
Here, a is the one-dimensional dark matter velocity dis- 
persion, set to (7 = Vq/^/2. We set the speed of the 
Sun around the Galactic center to be Vq = 220 km/s, 
for which the observational uncertainty is about 10% 
[53, 54]. The WIMP number density is = p-y^jm-^. 
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We assume! that the dark matter density is smooth and 
time-independent in the neighborhood of the Sun, and 
that = 0.3 GeV cm~^. Even if the dark matter were 
somewhat lumpy, the rcsuhs of the simulations will still 
be valid if is interpreted as the average density in the 
solar neighborhood [55]. 

Transforming to the heliocentric frame via a velocity 
transformation = v — v© , 

/s(x,v,)dWv, = //,(x,v, +Vo)d3xd3v,, (47) 

where the subscript s refers to quantities measured in 
the heliocentric frame. This distribution is anisotropic 
with respect to the plane of the solar system (the eclip- 
tic). The direction of the anisotropy with respect to the 
ecliptic depends on the phase of the Sun's orbit about 
the Galactic center. In order to avoid choosing a spe- 
cific direction for the anisotropy (i.e., to avoid choosing 
to start our simulations at a particular phase of the Sun's 
motion about the Galactic center), we angle- average this 
anisotropic distribution function to obtain an isotropic 
DF of the form 



fs{x,Vs) = j fs{y^,'Vs)dn 



(48) 



2(277)3/2 aVQVs 



-{vs-VQf/2a^ 



Using the angle-averaged DF is approximately valid for 
two reasons: (i) Scattering in the Sun is isotropic, so any 
bound WIMPs produced by elastic scattering will ini- 
tially be isotropically distributed, (ii) The time-averaged 
distribution function (averaged over the Sun's w 200 
Myr orbit about the Galactic center) only has a small 
anisotropic component [34], a consequence of the large 
(60°) inclination of the ecliptic pole with respect to the 
rotation axis of the Galaxy. If the flux at the Earth is 
dominated by particles whose lifetime in the solar sys- 
tem is greater than the period of the Sun's motion about 
the Galactic center, the use of time-averaged distribution 
function is appropriate. 

We use Liouville's theorem to find the halo DF for an 
arbitrary point in the Sun's potential well. Neglecting the 
gravitational potential of the planets and the extremely 
rare interactions between dark matter particles, each par- 
ticle's energy E is conserved: 



E = 



1 



(50) 
(51) 



where v is the speed of particle with respect to and in the 
gravitational sphere of influence of the Sun, and ^©(r) is 
the gravitational potential of the Sun ($0 = —GMq/t 
for r > Rq, where Rq represents the surface of the Sun). 
Thus, the DF within the Sun's potential well is 



f{'r,v) = fs{r,Vs{r,v)), 



(52) 
(53) 



An important consequence of this result is that the dis- 
tributi on functio n is identically zero for local velocities 
V < -\/— 2$0(r) = Vescir) below the escape velocity at 
that distance. 



B. Astrophysics Assumptions 

The Sun: The Sun is modeled as spherical and non- 
rotating. The gravitational potential and chemical com- 
position are described by the BS(OP) solar model [56]. 
We include the elements ^H, ^He, ^^c, ^'^m, ^^O, ^ONe, 
2^Mg, "^^Si, and ^^Fe in computing the elastic scattering 
rate. 

We treat the Sun with the "zero-temperature" approx- 
imation (i.e., the solar nuclei are at rest) since the kinetic 
energy of nuclei in the Sun is much less than the kinetic 
energy of dark matter particles. At the center of the Sun, 
the temperature is Tc ~ 10'' K, so the average kinetic en- 
ergy of a nucleus is of order 



1 keV. 



Ka =lkT, 



(54) 
(55) 



In the cooler outer layers of the Sun, the nuclei have even 
less kinetic energy. In contrast, the kinetic energy of dark 
matter particles in the Sun is of order 



100 GcV 



keV. 



(56) 
(57) 



The rms velocity of the nuclear species A is (i"^)^^^ = 
^2KA/mA ~ 500(mA/GcV)-^/2 j^j^^ g-i^ jQ^gj. ^j^^n the 
^ lO^ km/s speed of dark matter particles. Therefore, to 
good approximation, one can treat the baryonic species 
in the Sun as being at rest (i.e., having T = 0) 

The Solar System: The solar system is modeled as 
having only one planet, Jupiter, since Jupiter has the 
largest mass of any planet by a factor of 3.3 and there- 
fore dominates gravitational scattering. We place Jupiter 
on a circular orbit about the Sun, with a semi-major axis 
aiy_ = 5.203 AU, its current value, for the entire simula- 
tion, since its eccentricity is only e ~ 0.05 [57], and the 
fractional variation in its semi-major axis is < 10~^ over 
the lifetime of the solar system [58]. Jupiter is modeled 
as having constant mass density to simplify calculations 
of particle trajectories. This is not a realistic represen- 
tation of Jupiter's actual mass density but only a small 
fraction of particles scattered by Jupiter actually pene- 
trate the planet. WIMP-baryon interactions in Jupiter 
arc neglected since the optical depth of Jupiter is small 
enough that the probability of even one scatter occurring 
in the simulation is significantly less than unity. 

The orbit of Jupiter defines the reference plane for the 
simulation. The Earth's orbit is assumed to be coplanar 
with the reference plane, since the actual relative incli- 
nation of the two orbits is currently only 1.3°. 
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C. Initial Conditions 

The goal of this section is to compute the rate of elas- 
tic scattering of halo WIMPs by baryons in the Sun onto 
bound orbits, as a function of the energy and angular 
momentum of particles after the scatter. We also show 
how wc use this to choose the initial starting positions 
and velocities of the particles. There are two natural 
approaches to this problem: (i) Sample the dark matter 
flux through a shell a distance R > Rq from the cen- 
ter of the Sun, treating scatter in the Sun probabilisti- 
cally, and keeping only those particles which scatter onto 
Earth-crossing bound orbits, (ii) Calculate the scatter- 
ing probability in the Sun directly. The second approach 
is more efficient, and this is the one described below. 

The initial energy of a dark matter particle is 



(58) 

(59) 



where we have expressed the gravitational potential in 
terms of the local escape velocity from the Sun. The 
final energy of the dark matter particle is 



E' = E-Q 



/2 



(60) 
(61) 



where Q is the energy transfer between the dark matter 
particle and the nucleus during the collision. The energy 
transfer can be expressed in terms of the center-of-mass 
scattering angle 9 as (cf. Eq. A5) 



Qiv, cose) =2^v^ 
niA 



1 



COS( 



where 



Ma 



niA + m-^ 



(62) 



(63) 



is the reduced mass for a nuclear species of nucleon 
number A. The maximum energy transfer Qmax = 
2^.\v^/mA occurs if the dark matter particle is back- 
scattered, i.e., 9 = TT. Since we are interested in particles 
that scatter onto bound, Earth-crossing orbits,^ the in- 
teresting range of outgoing energy is 



' 2(0.5a®) 



< E' < 0, 



(64) 



where is the semi-major axis of the Earth's orbit, with 
the lower bound determined by the fact that the aphelion 
of a highly eccentric orbit is 2a. 

For a given incoming energy E, it may not be kine- 
matically possible to scatter into the full range of bound. 
Earth-crossing orbits. In particular, li E — Qmax = 
-E'min > 0' particle cannot scatter onto a bound orbit 
at all. Therefore, the lower bound on allowed outgoing 
energy is 



Emin = max 



GM, 



2(0.5a, 

while the upper bound remains 



, min(£; - Qmax{v), 0) , (65) 



EL 



0. 



Thus scattering rate of particles onto bound, 
crossing orbits is 



(66) 



Earth- 



Att r'^nA{r)v 



f{r,v) 



drdnrdvdQ ^ " ' dQ 

X e(J?0 - r)Q[-E'] 

X Q[E'-ELJ, (67) 

where we have imposed spherical symmetry on the Sun, 
Q.r is the solid angle in the Sun, /(.x, v) is the distribution 
function in Eq. (52), do-^i/dQ is the WIMP-nucleus cross 
section (Eq. Al), and 8(a;) is the step function. Since 



dE' = dQ, 
we can write Eq. (67) as 



- = 477 2_^r nA{r)v''-^f{r,v) 



(68) 



drdQ.rdvdE' 



dQ 

X e(i?o - r)Q[-E'] 

X e[E'-ELi^], (69) 



By sampling this distribution, we find the initial energy, 
speed, and scattering position vector of the WIMPs. 

The outgoing energy is distributed uniformly unless 
there is kinematic suppression. The kinematic suppres- 
sion is most pronounced for large WIMP masses and very 
negative energies because, in order for a particle to scat- 
ter onto a bound orbit. 



Vs < 2— Vesc{r) 

rriy — niA 



(70) 



In principle, particles scattered to bound orbits with a < a0/2 
could later evolve onto Earth-crossing orbits. However, the 
torque from Jupiter is never high enough to pull a particle with 
a < 09/2 onto an Earth-crossing orbit unless — ci)/a ^ 

10~^. Moreover, each additional scatter in the Sun reduces the 
energy of the orbit in the limit of a cold Sun, so the semi-major 
axis may only shrink. 



where is the particle velocity at infinity. Heavy dark 
matter particles can only scatter onto bound orbits if 
their velocities at infinity are only a small fraction of 
the escape velocity from the Sun a distance r from the 
Sun. For energies for which the kinematic suppression 
is minimal, we express the uniformity of diV0/d£" in 
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terms of the semi-major axis. Since E' = —GMQ/2a for 
particles on elliptical orbits, dN^^/da cx a^^, or 



dloga 



-1. 



(71) 



Therefore, most particles scatter onto relatively small or- 
bits. 

The angular momentum of each scattered particle is 
in the range J € [0,rv'], where r is the radius from the 
center of the Sun at which the particle scatters. To de- 
termine the distribution of magnitudes and directions for 
the angular momentum, we assume that the direction of 
the final velocity v' is isotropically distributed with re- 
spect to the position vector r. If we specify 9^ to be the 
colatitude of the velocity vector with respect to the po- 
sition vector, and the magnitude of the angular momen- 
tum is J = rv' sin 9y^ then the distribution in angular 
momentum at fixed r, v' has the form 



dNm oc d cos 



dj2 



2r2w'2 ^1 - J2/(r2i;'2) 



(72) 



The effect of kinematic suppression due to a large WIMP 
mass is that the particles that do scatter onto bound or- 
bits can only do so close to the center of the Sun where 
Vesc is highest. This reduces the maximum angular mo- 
mentum of the outgoing particle, and so eccentricity is 
an increasing function of WIMP mass. 

By sampling Eq. (69), Eq. (72), and the azimuth of v' 
with respect to the position vector, we fully specify the 
6-dimensional initial conditions of the WIMPs, sampled 
to the same density as they would actually scatter in the 
Sun. 



D. Simulation Specifics 

The goals of the simulations are to predict the direct 
and indirect detection signals from particles bound to the 
solar system (relative to the signal from unbound parti- 
cles) as a function of m-^ and the elastic scattering cross 
section. The simulate orbits for a variety of WIMP pa- 
rameters and then interpolate the results for other values 
of those parameters. 

We ran four sets of simulations with different choices of 
and . The first simulation, called "DAM A", used 
= 60 atomic mass units (AMUs) and a^^ — 10"''^ 



cm . These parameters lie in the DAM A annual mod- 
ulation region [59, 60]. A second simulation, called 
"CDMS" , used the same WIMP mass as in the DAMA 
simulation but a cross section two orders of magnitude 
lower, cTp^ = 10~^^ cm2, below the minimum of the 
CDMS 2006 exclusion curve (Fig. 4). Two more sim- 
ulations were chosen to have — 10"'*^ cm2 but with 
larger WIMP masses. The "Medium Mass" simulation 
uses = 150 AMU, and the "Large Mass" simulation 
uses rriy = 500 AMU, selected to explore the dependence 
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FIG. 4: Points in the err 



parameter space used for the 



solar capture simulations, plotted along with exclusion curves 
from recent experiments. This plot was generated with 
http: / /dendera.berkeley.edu / plotter / entryform.html . 



of the simulations on WIMP mass. The choices for 
and <7p^ are plotted in Fig. 4 in addition to some re- 
cent direct detection results. The details on the initial 
conditions of the simulations are summarized in Table I. 

Since integrating the orbits of particles in the solar sys- 
tem is computationally expensive, it is more important 
to integrate just enough orbits to determine the approx- 
imate size of the bound DF relative to the unbound dis- 
tribution, and to get a sense of which effects matter the 
most, than it is to get small error bars on the DF. This 
principle guides our choices in the sizes of the ensembles 
of particles. 

The number of particles simulated Np in each of the 
solar capture simulations is given in Table I. We follow 
particles with semi-major axes slightly below the Earth- 
crossing threshold so that if the semi-major axis increases 
modestly during the simulation, the contribution to the 
Earth-crossing flux is included. Fewer particles were sim- 
ulated in the runs with lower cross section because the 
lifetimes were far longer than in the DAMA run. 

We use the flowchart in Fig. 5 to sketch the simu- 
lation algorithm. There are six things that need to be 
set in order to run the simulations: starting conditions; 
the radius Vc at which the heliocentric-barycentric coor- 
dinate change needs to be made; methods for initializing 
h and potentially changing h at later times; conditions 
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TABLE I: Solar capture simulations 



Name 


[AMU] 




cm^] 


Np [a > 0.48 AU] 


DAMA 


60 


10 


-41 


1078586 


CDMS 


60 


10 


-43 


145223 


Medium Mass 


150 


10 


-43 


144145 


Large Mass 


500 


10 


-43 


148173 



for passing through and scattering in the Sun; the size of 
the bubble about Jupiter, Iq^, and the accuracy criterion 
\AE/E\: and conditions for terminating the simulation. 
Following the flowchart, we discuss each point in turn. 

Starting Conditions We sample the distribution of 
WIMPs initially scattered into the solar system accord- 
ing to Eqs. (69) and (72). Once we have determined 
the initial position and velocity of each WIMP, we trace 
the WIMPs back to perihelion and start the integration 
there. We follow all particles after the initial scatter, us- 
ing the map to evolve the particles to the Sun bubble 
wall (0.1 AU) using map described in Section II CI. In 
order to account for the fact that particles may experi- 
ence a second scatter before leaving the Sun for the first 
time, we perform a rescattering Monte Carlo when we 
construct the DFs. 

Once the particles have reached the bubble boundary, 
we initialize the adaptive time step symplectic integra- 
tor (Section II A), setting h = 10~^i?g^ yr and integrat- 
ing the equations of motion in heliocentric coordinates. 
With this choice of initial h, a particle with initial semi- 
major axis = 1 AU will be integrated with 4.7 x 10^ 
steps/orbit, while a particle with a = 100 AU will be 
integrated with 4.7 x lO*' steps/orbit. We choose such 
a small h to minimize errors near perihelion, which is 
the point in the orbit at which errors are largest (Sec- 
tion II B). If the semi-major axis exceeds rc/2, it may 
be necessary to change to barycentric coordinates before 
the particle reaches aphelion for the first time. 

Coordinate Change For the weak scattering simula- 
tions, we set e = 0.1 (Eq. 42), thus setting the transition 
radius between the heliocentric and barycentric coordi- 
nated regimes to rc = 53 AU. This is large enough that 
only a small percentage of particles routinely cross the 
transition radius, but small enough that the heliocentric 
potential does not have too large a contribution from the 
indirect potential. 

Setting h: After the particles reach their first aphelion, 
h is reset according the values listed in Table II. These 
values of h are chosen such that both errors at perihelion 
{\/^E/E\ < IQ-^) and near Jupiter {\/:^E / E\ < 10"^) 
are small. The combination of the values of h and the 
Jupiter bubble radius Iqy. (see below) were determined 
empirically. We used slightly smaller values of h for some 
semi- major axes in the CDMS, Medium Mass, and Large 
Mass runs compared to the DAMA run in order to check 



TABLE II: The fictitious time step /i as a function of semi- 
major axis a for the DAMA simulation and the simulations 
with a^^ = 10"*^ cm^ (CDMS, Medium Mass, Large Mass). 



a [AU] 


DAMA 

h [RJ;' yr] 


erf = 10-« cm2 

h [R^,' yr] 


a < 0.75 


1 


X 10"" 


1 


X 10"" 


0.75 < a < 1.1 


7 


X 10"^ 


7 


X 10"® 


1.1 < a < 1.6 


6 


X 10"''' 


6 


X 10"® 


1.6 < a < 2.2 


5 


X 10"^ 


2 


X 10"® 


2.2 < a < 3.5 


4 


X 10"^ 


2 


X 10"® 


3.5 < a < 6.2 


3 


X 10"^ 


1.5 


X 10"® 


6.2 < a < 11 


2 


X 10"^ 


1 


X 10"® 


11 < a < 13 


9 


X lO"'' 


2 


X 10"® 


13 < a < 22 


2 


X 10"^ 


2 


X 10"® 


22 < a < 30 


2 


X 10"'^' 


2 


X IQ-^' 


30 < a < 45 


1 


X lO"'' 


1 


X 10"® 


45 < a < 120 


6 


X 10"^ 


6 


X 10"^ 


120 < a < 200 


4 


X 10"'^ 


4 


X 10"'' 


200 < a < 500 


3 


X 10"'^ 


3 


X 10"'' 


a > 500 or unbound 


2 


X 10"^ 


2 


X 10"'' 



that the behavior of long-lived WIMPs was not an arti- 
fact of the choice of parameters. 

A particle's energy (and hence, semi-major axis) may 
change throughout the simulation. If the energy changes 
by 20% or more from when the particle enters the Jupiter 
bubble to when it exits, the particle is flagged to have h 
adjusted at the next aphelion. We do not readjust h af- 
ter every aphelion, or after each time the particle passes 
through the bubble, because very frequent changes in h 
can induce growing numerical errors in the Jacobi con- 
stant. We impose any changes in h at aphelion instead 
of the bubble boundary, since we have determined exper- 
imentally that aphelion is the optimal point at which to 
change h. 

The Sun Bubble When a particle first crosses into the 
bubble about the Sun, we calculate its perihehon. If the 
perihelion is smaller than 2Rq, we proceed to map its 
current position and velocity to its position and velocity 
as it exits the bubble according to the prescription of 
Section II CI. If the perihelion lies within the Sun, we 
employ a Monte Carlo simulation of scattering in the 
Sun (Appendix B). The vast majority of the time, the 
particle does not scatter, and we simply use the map to 
move the particle to the edge of the bubble. If the particle 
rescatters onto an orbit with semi-major axis a < 0.3 AU, 
we terminate the integration. If the particle rescatters 
onto an orbit with o > 0.3 AU, we evolve the new orbit 
to the edge of the bubble, and then resume the adaptive 
time step symplectic integration. 

The Jupiter Bubble For the DAMA, CDMS, and 
Medium Mass simulations, we set the Jupiter bubble 
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FIG. 5: Flowchart for the simulation algorithm for the solar capture experiments. 



boundary to be ii,. = 1.7 AU, and the accuracy crite- 
rion to be \IS.Ef/Ef\ < 10^^. We adjusted this value 
for some particles in order to speed up the integration 
in cases where particles had generically slightly smaller 
initial \AEi/Ei\ than \AE f / Ef\, and took a longer time 
with \ AEf/Ef \ = 10~^ to reach a sufficiently flat plateau 
in \A.Ef/Ef \ than with a slightly larger accuracy crite- 
rion. Wc set Iq^ = 3.4 AU past 500 Myr for the Medium 
Mass simulation to determine if there were any effects of 
a larger Zq^ on the orbits. There were no effects on the DF 
resulting from this change. For the Large Mass simula- 
tion, we experimented with a lower value of the accuracy 



criterion ( \AEflEf \ < 2 x 10"'^ for the first 500 Myr, 
\AEf/Ef \ < 3x10"^ later) and a larger bubble, \ = 2.1 
AU. The only effect the larger accuracy criterion had was 
to raise slightly the maximum energy error per orbit. 

Stopping Conditions There are three reasons for ter- 
minating an integration. If the particle crosses outward 
through the shell r — 5000 AU, the integration stops. 
Particles crossing this shell will rarely cross the Earth's 
orbit again. Secondly, if the particle rescatters in the Sun 
onto an orbit with a < 0.3 AU, we halt the integration 
since the particle will soon thermalize in the Sun and 
will not cross the Earth's orbit. Thirdly, we stop the in- 
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tegration if the particle survives for a time = 4.5 Gyr, 
approximately the lifetime of the solar system. 



E. Computing 

Simulations were performed on three Linux beowulf 
computing clusters at Princeton University. Each set of 
simulations required 10^ CPU cycles on 3 GHz dual core 
processors. 



IV. DISTRIBUTION FUNCTIONS 

Before presenting the results of the simulations, we de- 
fine some terms that will be used frequently in this sec- 
tion. The "heliocentric frame" describes an inertial frame 
moving with the Sun. "Heliocentric speeds" will refer 
to speeds relative to the Sun, measured at the Earth as- 
suming the Earth has zero mass. The "geocentric frame" 
refers to an inertial frame moving with the Earth. Unless 
otherwise noted, geocentric WIMP speeds are those out- 
side the potential well of the Earth. The "free space" dis- 
tribution fimction, Eq. (49), is the anglc-avcragcd halo 
distribution function in an inertial frame moving with 
the Sun, outside of the gravitational sphere of influence 
of the Sim. The "unbound" distribution function refers 
to the Liouville transformation of the free space distri- 
bution function to the position of the Earth (Eq. 52), 
including the effects of the gravitational field of the Sun 
but not the Earth. 

In Fig. 6, we present the one- dimensional geocen- 
tric DFs divided by the halo WIMP number density 
(defined in Section HI A) for each simulation. These 
DFs have already been integrated over angles, and are 
normalized such that the bound dark matter density 
%,6oMnd = / dvv^fiv), where f{v) = J dl^/(v). We 
plot the DFs in Fig. 6(a) on a logarithmic scale in or- 
der to highlight the structures, while we plot the CDMS 
simulation (Table I) DF on a linear scale in Fig. 6(b) 
to compare the simulation results with theoretical DFs. 
The DFs are based on (1 — 4) x 10^ passages of parti- 
cles within a height \z\ < Zc = lOi?© of the Earth's orbit, 
and estimated using the technique described in Appendix 
C. The DFs are insensitive to Zc, at least in the range 
1 i$ -^c ^ 25i?0. Error bars are based on 500 bootstrap 
resamplings of the initial scattered particle distributions 
for each simulation. 

The most striking feature of the DFs is the smallness 
with respect to the DF of halo WIMPs unbound to the 
solar system. This is in stark contrast to the prediction 
of Damour and Krauss [37] . In order to show why this is 
the case, we examine the simulations in detail. In partic- 
ular, we (i) identify the features in the DF with specific 
types of orbits, (ii) find the lifetime distribution of such 
orbits, and (iii) determine what sets the lifetime of those 
orbits (e.g., ejection vs. rescattering and thermalization 
in the Sun). With these data, we may also determine 



how the DF varies with WIMP mass and cross section, 
and estimate the maximum DF consistent with limits on 
the spin-independent cross section. 

The DFs from the four simulations show similar mor- 
phologies, although the normalization of the features dif- 
fers. The most prominent feature in all four DFs in Fig. 
6 is the "high plateau" between 27 < f < 48 km s^^. 
In order to identify which orbits contribute to this 
plateau, it is useful to examine the two-dimensional DF. 
In Fig. 7, we show the two-dimensional DF f{v,cos9) = 
J A(l)f{v, cos e,cj)) for the CDMS simulation (Table I) in 
both (a) heliocentric and (b) geocentric coordinates. The 
angle between the velocity vector and the direction of 
the Earth's motion is 9, while (j) is an azimuthal angle, 
with (j) — Q corresponding to the direction of the north 
ecliptic pole if 9 = it/2. The DFs are plotted on a log- 
arithmic scale to highlight structure. We only show the 
CDMS simulation results in this figure since the phase 
space structure of the DF is virtually the same in all 
simulations. 

From Fig. 7(b), we identify the short arc between 
27 < w < 50 km s~^ below cos9 < —0.5 with the high 
plateau, although there is a small contribution from the 
other, longer arc. We find that the short arc in the geo- 
centric DF corresponds to the trumpet-shaped feature in 
the heliocentric DF below w < 38 km s~^. For bound 
orbits, the heliocentric speed at r = 1 AU is 



v{a) 



1/2 



(73) 



where is the orbital speed of the Earth. The helio- 
centric speed w = 38 km s~^ corresponds to the lowest 
Jupiter-crossing orbit, so the trumpet feature in the two- 
dimensional heliocentric DF corresponds to orbits that 
do not cross Jupiter's orbit. 

The trumpet shape of the heliocentric DF (and the 
narrow band in the geocentric DF) in Fig. 7 can be sim- 
ply explained. In Fig. 8, we calculate the energy and 
angular momentum for each point in velocity space. The 
black region of velocity space represents unbound orbits. 
All points for which orbits are bound and have perihe- 
lia inside the Sun are marked in dark grey, while orbits 
that are bound and cross Jupiter's orbit are marked in 
light grey. The white regions correspond to bound or- 
bits that neither enter the Sun nor cross Jupiter's path. 
If we were to integrate the ^-slices shown if Fig. 8, we 
would find that the region in Fig. 8 corresponding to 
Sun-penetrating orbits that do not cross Jupiter exactly 
matches the parts of phase space we identified with the 
high plateau. 

Fig. 8 was computed for a system without planets. 
Once Jupiter is added, another type of orbit that may 
exist is a bound orbit for which is fixed but J is not. 
An example of this type of orbit is a Kozai cycle. In 
this case, = a^v cos 9 in the heliocentric frame. In the 
special case that (f) = t: /2, J = J^. Therefore, the parts of 
{v, cos9) phase space in the = 7r/2 plane corresponding 
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FIG. 6: Geocentric distribution functions from the simulations, (a) Results from all simulations, (b) The CDMS distribution 
function relative to theoretical distribution functions for unbound WIMPs. 
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FIG. 7; Distribution functions divided by in the v~cos8 plane (integrated over (p) for both (a) heliocentric and (b) geocentric 
frames. These DFs come from the CDMS simulation, and the units are (km s~^)~^ 



to Sun-penetrating orbits also cover orbits witli fixed 
by the initial scatter in the Sun for other values of (/). 
Thus, the high plateau in the 1-dimensional geocentric 
DF is built up by WIMPs with a < and small 

but not necessarily small J. 

The second feature of the distribution functions in Fig. 
6 is the "low plateau." This is the relatively flat part of 
the distribution function that extends from « 10 km s^^ 
to « 70 km s^^. This feature corresponds in the long arc 
in the two-dimension DF in Fig. 7 and the stripe between 
38 < w < 42 km s~^ in the heliocentric DF. From Fig. 



8, we identify this feature with bound, Jupiter-crossing 
orbits. Small gaps exist in the heliocentric DF with 
w > 40 km s^^ and cos0 < and 38 < < 40 km s^^ 
and cos 9 > because these regions of phase space are 
inaccessible to WIMPs initially scattered in the Sun in 
the restricted three-body problem. This translates to a 
truncation of the low plateau at the extrema in geocentric 
speed. 

The third common set of features in the one- 
dimensional geocentric distribution function are spikes 
in the low plateau. These spikes result from the long- 
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lifetime tail of Jupiter-crossing or nearly Jupiter-crossing 
particles that spend significant time near mean-motion 
resonances or on Kozai cycles. The minimum semi- 
major axis for these spikes corresponds to the 3:1 reso- 
nance, a « 2.5 AU. Long lifetime tails due to "resonance- 
sticking" orbits have also been observed in simulations of 
comets [61, 62]. The error bars on the spikes are large 
due to the small numbers of long-lived resonance-sticking 
particles in each simulation. There is also some Poisson 
noise in the height of the spikes between simulations due 
to the small number long-lived WIMPs in each simula- 
tion. 

Next, we show the lifetime distribution of bound 
WIMPs and demonstrate which mechanisms (thcrmal- 
ization or ejection) terminates the contribution of orbits 
to the DP. In Fig. 9, we present the lifetime distribu- 
tions for all WIMPs in the DAMA, CDMS, Medium Mass 
and Large Mass simulations. There are several notable 
features in this plot. First, and most striking, many of 
the bound particles survive for very long times — up to 
10® — lO'* yr. However, in none of the simulations is there 
a large population of particles that survive for times of 
order the age of the solar system, although there is a 
small population that does (the notch at 4.5 Gyr in Fig. 
9). Secondly, the lifetime distribution functions of the 
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FIG. 9: Particle lifetime distributions for the DAMA (solid), 
CDMS (dot-dashes), Medium Mass (short dashes), and Large 
Mass [long dashes) simulations. 
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CDMS, Medium Mass, and Large Mass runs are very 
similar. However, these lifetime distributions are quite 
different from that of the DAMA simulation. The differ- 
ences in the lifetime distributions are due almost entirely 
to the elastic scattering cross section, at least for the 
range of WIMP masses we consider. 

In order to both explain these features in the lifetime 
and phase space distribution functions, it is useful to ex- 
amine the lifetime distributions as a function of the initial 
semi- major axis a^, as shown in Fig. 10. 

The largest feature in Fig. 9 is the strong peak at 
ti ~ 10^ yr for the DAMA simulation and U ~ 10^ yr 
for the simulations with a^^ — 10^^'^ cm^, which we 
call the "rescattering peak." It encompasses the major- 
ity of particles in each simulation. This feature results 
from WIMPs that rescatter in the Sun before they are 
ejected from the solar system by Jupiter or precess onto 
orbits that do not intersect the Sun. This rescattering 
peak is offset between DAMA and the other simulations 
because the lifetime is inversely proportional to the scat- 
tering probability in the Sun, ti cx (J~^. 

There is one important difference in the composition 
of the rescattering peak between the DAMA and other 
simulations. In Fig. 10, we show that particles on 
Jupiter-crossing orbits exhibit a rescattering feature in 
the DAMA simulation but not in the other simulations. 
Indeed, about 23% of Jupiter-crossing particles in the 
DAMA simulation are rescattered in the Sun, while < 2% 
are rescattered in the other simulations. This is because 
the timescale on which Jupiter can perturb the perihelia 
of Jupiter-crossing orbits out of the Sun is significantly 
shorter than rescattering timescales for the a^^ = 10"^'^ 
cm^ simulations, but the two timescales become closer at 
higher cross sections. 

Another feature occurs at ti ^ 10^ yr, which we call the 
"ejection peak." This feature occurs at the same time for 
each simulation, and from Fig. 10 we see that this arises 
from Jupiter-crossing orbits. The median time at which 
this feature occurs is independent of (7^^ since ejection, 
not rescattering, sets the lifetime of these WIMPs. The 
slope of the Jupiter-crossing lifetime distribiition changes 
near ^ 10 Myr for all simulations. WIMPs that have 
ti > 10 Myr are resonance-sticking, and their lifetime 
distribution goes as N{t) oc where a is slightly less 
than one. 

A third feature, called the "quasi-Kozai peak," is cen- 
tered at ti ~ 10^ yr in the DAMA simulation, and 
ti ~ 10* yr in the other simulations. The feature is seen 
in the 1.5 AU < a, < 2 AU and 2 AU < a, < 
bins of Fig. 10. The WIMPs in the quasi-Kozai peak 
are not on true Kozai cycles because of significant in- 
teractions with mean-motion resonances. In the simu- 
lations, particles in the quasi-Kozai peak are observed 
to alternate between rescattering peak-type orbits with 
perihelion well inside the Sun, and orbits that look like 
Kozai cycles. Both the semi-major axis and Jz vary with 
time; neither is conserved although the combination giv- 
ing the Jacobi constant Cj is fixed (Eq. 20). There are 



some orbits at the low end of the semi-major axis range 
1.5 AU < a < for which a and Jz are conserved 

and the Kozai description is accurate. 

The median lifetime of WIMPs on quasi-Kozai cycles is 
well-described by ti/P^ « 300/r, where is the WIMP 
orbital period and r is the optical depth through the 
center of the Sun. This implies that particles are eventu- 
ally removed from Earth-crossing orbits by rescattering 
in the Sun. The height of the rescattering peak relative 
to the quasi-Kozai peak is greater in the DAMA simula- 
tion than the other simulations because the optical depth 
in the Sun is large enough that particles originating deep 
within the Sun rescatter before the torque from Jupiter 
can pull the perihelion towards the surface of the Sun. 

The fourth feature is not obvious in Fig. 9, but is once 
the lifetime distribution is displayed on logarithmic scales 
in Fig. 10. This feature is the "Kozai peak." This peak 
is located at about ti ~ 10* yr for the DAMA simula- 
tion, and near ~ = 4.5 Gyr for the other simula- 
tions. This feature results from particles whose orbital 
evolution can be described by Kozai cycles (a, Jz con- 
served), of which we see both circulating and librating 
populations. For the CDMS, Medium Mass, and Large 
Mass simulations, the peak is at tQ because that is the 
time at which we terminate the simulations. Particles 
on these orbits have a; < 1.5 AU, and originate in the 
outer r > 0.5i?© in the Sun. They constitute only a 
small fraction (~ 0.1%) of all orbits with at < 1.5 AU, 
but dominate the lifetime distribution of particles with 
lifetimes > 10* yr. The median lifetime of particles on 
Kozai cycles depends on the WIMP-nucleon cross section 
in the form ti/P^ « 10^/r, where again r is the optical 
depth through the very center of the Sun (t ~ 10~^ for 
DAMA, T ~ 10~^ for the other simulations). 

Now that we have identified the types and lifetimes of 
bound WIMP orbits, we see how these come together to 
build up the phase space DF as a function of WIMP mass 
and cross section. 



A. The Distribution Function as a Function of cr^ 

We have identified (i) which range of initial semi-major 
ax;is corresponds to the features in the geocentric DFs 
in Fig. 6, and (ii) described the lifetime distributions 
of WIMPs. Next, we describe the composition (not just 
in terms of the semi-major axis Oj but by the type of 
orbit) and height of the DFs as a function of ct^^. This 
is most easily illustrated with the time-evolution of the 
DFs, which we show in Fig. 11 for the (a) DAMA and (b) 
CDMS simulations. We focus on these two simulations 
since the salient results of the Medium Mass and Large 
Mass simulations closely resemble those of the CDMS 
run. In each plot, we show distribution functions as a 
function of time since the birth of the solar system, for 
t = 10^ yr, t = 10"^ yr, t = 10* yr, t = 10^ yr, and 
t = tQ= 4.5 Gyr. 

The low plateau, composed of Jupiter-crossing parti- 
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FIG. 10: Lifetime distributions as a function of initial semi-major axis. 



cles, has reached equilibrium in both ~ 10^ year for both 
simulations. The only growth in the low plateau after 
10 Myr comes from particles on resonance-sticking or- 
bits that pump up the spikes. The time evolution of the 
low plateau (but not its final equilibrium height) is in- 
dependent of cross section over two orders of magnitude 
in WIMP-baryon cross section because the equilibrium 
timescale is essentially the ejection timescale. The height 
of the low plateau is proportional to the rate at which 
particles are initially scattered onto Jupiter-crossing or- 
bits, Nq^. Since the scattering rate is proportional to 
the cross section, the height of the low plateau is pro- 
portional to the spin-independent cross section, at least 
in these simulations. One would expect that the plateau 
height would grow less rapidly with tTp ^ if the lifetimes of 
Jupiter-crossing WIMPs were dominated by rescattering 
in the Sun, not ejection from the solar system. 

The absolute height of the spikes is similarly related to 



Nif_ and the relative ejection and rescattering timescales; 
the spikes in the DAMA simulation are more prominent 
than in the CDMS simulation because N% is two orders 
of magnitude larger. The time-evolution of the spikes 
can be explained by the following. The lifetime dis- 
tribution of spike WIMPs falls as N{t) cx where 
a is slightly less than one. The rate at which WIMPs 
cross the Earth's orbit is Nc{t) = const if the long-lived 
WIMPs are resonance-sticking. Therefore, the total con- 
tribution of the spike WIMPs to the DF beyond time t 
goes as 

fsp^kei>t) cc j ° N{t')N,{t')At' (74) 

oc t;^-"-ti-". (75) 

This argument is only strictly true if the types of spike 
orbits is independent of the lifetime distribution, which is 
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FIG. 11; Growth of the distribution functions as a function 
of time for the (a) DAMA and (b) CDMS simulations. 



uncertain due to the small number statistics of the spike 
WIMPs. However, in Fig. 11, some of the spikes either 
grow linearly with time or do not grow at all for some 
stretches of time. This phenomenon is due to the small 
numbers of long-lived resonance-sticking particles. For 
an individual WIMP, f{v) (x t ii t < ti and f{v) is fixed 
for t > ti. 

The time evolution of the high plateau is different be- 
tween the DAMA and CDMS simulations. In the DAMA 
simulation, by i = 1 Myr, nearly all of the rescattering 
peak WIMPs have rescattered and thermalized in the 
Sun, as have a significant fraction of the quasi-Kozai 
WIMPs. At this point, the high plateau in the range 
27 km s~^ < u < 45 km s"""^ is built up by roughly 
similar contributions from rescattering and Kozai peak 
WIMPs. The contribution of Kozai WIMPs relative to 
rescattering peak WIMPs at a particular time can be es- 
timated using 

~ (3 (-^) , (76) 

Jrescatt \'^med / 



where f Kozai is the DF due to Kozai cycling WIMPs, 
frescatt that of rcscattcring peak WIMPs, /3 is the frac- 
tion of WIMPs with a < 1.5 AU on Kozai cycles, and 
t„ied is the median lifetime of rescattering peak WIMPs. 
f3 « 10~^ for WIMPs experiencing Kozai cycles, and 

tmed ~ 10^ yr, so f Kozai/ frescatt ^ I at t ^ 1 Myr. 

This assumes that the increase in the DF for a Kozai cy- 
cling WIMP as a function of time is similar to that of a 
WIMP on a rescattering peak- type orbit. The feature in 
the DF between 45 km s~^ < w < 50 km s~^ is due to 
quasi-Kozai WIMPs. 

The high plateau grows substantially between t = 1 
Myr and t = 100 Myr, although not strictly linearly be- 
cause scatters in the Sun remove Kozai and quasi-Kozai 
WIMPs from Earth-crossing orbits. The error bars on 
the DF increase with time as the ever-decreasing num- 
ber of Earth-crossing WIMPs (Fig. 10) build up the DF. 
After 100 Myr, the high plateau grows very slowly until 
it reaches equilibrium hy t = \ Gyr; in our simulation 
of 8 X 10^ WIMPs with orbits interior to Jupiter's orbit, 
only one WIMP has a lifetime of 1 Gyr. 

Even though we simulate ~ 10^ WIMPs on Kozai cy- 
cles, we are clearly undersampling those with ti > 10^ 
yr. To estimate how much larger the DF could be, we 
note that the lifetime distribution of Kozai WIMPs with 
ti > 100 Myr is well fit by N{t) oc t^^. If we assume that 
the rate at which the long-lived WIMPs contribute to 
the DF as a function of time is the same as for the Kozai 
WIMPs we simulate, then Nc{t) — const. Therefore, 
according to Eq. (74), the part of the DF built after 
time t is fKozai{v,> t) oc t^^ . For the high plateau, 
fKozai{v,t > 10* yr) > fKozai{v,t > 10^ yr), so we be- 
lieve that we have not underestimated the high plateau. 

A major consequence of this equilibrium distribution 
is that the high plateau of the DF f{v)/n^ is fixed es- 
sentially fixed above a certain cross section. Since the 
lifetime of Kozai orbits U cx {(Jp^)^^, we find that the 
high plateau is is fixed for (7^^ > lO"**^ cm^. 

The time evolution of the distribution function is a 
bit different in the simulations for which a^^ = lO^^'^ 
cm^. At t = 1 Myr, the high plateau is dominated by 
rescattering peak WIMPs, which have a median lifetime 
tmed ~ 10^ yr. Between t = 1 Myr and t = 100 Myr, the 
growth in the high plateau is mostly due to the long- 
lifetime tail of the rescattering peak WIMPs and the 
quasi-Kozai WIMPs. This is because f Kozai/ frescatt ~ 1 
only for t ^ 10* yr. For t > 100 Myr, the high plateau 
is dominated by Kozai WIMPs. To determine if we suffi- 
ciently sampled the Kozai population, we compared the 
DF derived from the DAMA simulation when the in- 
tegrated optical depths were equivalent to those in the 
CDMS simulation (if effect, comparing the DAMA simu- 
lation a,t t — 10^ yr with the CDMS simulation a.t t — 10^ 
yr). We found the DFs to be consistent with each other. 

Unlike in the DAMA simulation, a number of particles 
have lifetimes longer than the age of the solar system 
(~ 100 out of ^ 10"''). One consequence of this is that the 
DFs should be somewhat smaller than the DAMA distri- 
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bution function, since the high plateau of the DAMA sim- 
ulation has reached equilibrium by the present but the 
a^^ — 10"^'^ cm^ distribution functions are still grow- 
ing. In fact, we find that the high plateau is about a 
factor of three smaller for the CDMS simulation than for 
the DAMA simulation. As decreases, so should the 
height of the high plateau. For < 10~^^ cm^, the 
high plateau should be dominated by rescattering peak 
orbits. 

In summary, we find that while the DF for cr^ = 10^^^ 
cm^ is dominated by Kozai WIMPs, there is some con- 
tribution from long-lived Jupiter-crossing WIMPs (al- 
though the error bars are large due to small number 
statistics). As the cross section decreases, the Jupiter- 
crossing component of the number density also decreases, 
and the Kozai and quasi-Kozai contributions dominate. 
However, the Kozai WIMPs fail to reach equilibrium, so 
the overall DF goes down as a function of decreasing 
cross section. Below ~ 10~^^ cm^, we expect the DF to 
be dominated by the rescattering peak WIMPs. 
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FIG. 12: Percentages of particles in each initial perihelion bin. 
Poisson errors smaller than points. 



There is little variation in the morphology of the life- 
time distributions for the three simulations with a^^ = 
10"''^ cm^. The shape of the lifetime distribution appears 
to be determined almost solely by the elastic scattering 
cross section, not the particle mass, at least in the range 
of masses considered in these simulations. It is possible 
that these distributions (in lifetime and density compo- 
sition) for a very high or very low mass WIMP would 
perhaps look different from those in Fig. 10. 

However, the DFs in Fig. 6 did show some variation 
with WIMP mass. There are three effects that might 
induce a mass dependence on the DF. (i) The mass can 
affect the initial energy and angular momentum distribu- 
tion of bound WIMPs. As discussed in Section III C, it is 
increasingly difficult to scatter halo WIMPs onto bound 
orbits as the WIMP mass increases. The maximum en- 
ergy transfer Qmax approaches an asymptote for large 
WIMP masses, but the unbound WIMP energy increases 
since energy E cx m^. Thus, the minimum scattered 
particle energy E' — E ~ Qmin increases for fixed initial 
speed but increasing WIMP mass. However, this is not 
expected to be a major effect for the range of masses used 
in the simulations. 

The angular momentum distribution is also affected by 
the WIMP mass, as parametrized by the initial particle 
perihelion in Fig. 12. As discussed in Section IIIC, the 
maximum angular momentum decreases with increasing 
since high mass particles scattering onto bound or- 
bits must do so at smaller distances from the center of 
the Sun. Thus, the Medium Mass and Large Mass simu- 
lations have a deficit of large perihelion particles relative 
to the CDMS simulation. Since Kozai WIMPs originate 
in the outskirts of the Sun, this suggests that there will 
be fewer particles on Kozai cycles as the WIMP mass 



increases. 

(ii) The particle mass affects the rescattering prob- 
ability in the Sun. In Eq. (B8), we show that the 
scattering probability along a path I is proportional to 
dr/dZ oc (^1 — e^'^'""^/'^^) , which is a mildly increas- 
ing function of WIMP mass (since Qmax is mass- 
dependent, Eq. 62). The optical depth for the Large 
Mass simulation (m^ — 500 AMU) for a given path is 
about 15% higher than for rriy^ = 60 AMU. However, 
while high mass WIMPs have a higher scattering prob- 
ability than low mass WIMPs, they also rescatter far 
more often onto Earth-crossing orbits. Therefore, it is 
not clear from the outset whether high mass WIMPs 
will have longer or shorter lifetimes relative to low mass 
WIMPs. 

(iii) The WIMP mass also affects the overall amplitude 
of the final bound dark matter DF because the WIMP 
mass determines the scattering rate of halo particles onto 
bound. Earth-crossing orbits. For high mass WIMPs, the 
total capture rate of halo WIMPs in the Sun is [e.g., 63] 

Ntot/n^ m~^, m^:^mA- (77) 

The function Ntot/n^ is plotted in Fig. 13(a) for the 
capture rate due to all species in the Sun {solid red) 
and for scattering only on hydrogen {blue dots; calcu- 
lated in the limit of a cold Sun). The capture rate 
of particles onto Earth- crossing orbits is shown in Fig. 
13(b). Note that the capture rate N^/uy, onto Earth- 
crossing orbits is an increasing function of WIMP mass 
until about « TeV (« 100 GeV in the case of only 
hydrogen scattering). This is because low mass WIMPs 
may be scattered onto very small orbits (whose aphelia 
may be within the Sun), which are kinematically sup- 
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pressed for higher mass WIMPs. Even though the total 
WIMP capture rate decreases for higher WIMP mass, 
those WIMPs that are captured are increasingly prefer- 
entially scattered onto Earth-crossing orbits. The func- 
tion N^/n^ turns over when most captured particles are 
on Earth- crossing orbits, and then the function follows 
the familiar N^/n^ oc 'rn~^ . 

The consequence of these scattering rates of halo 
WIMPs in the Sun is that, if the DFs were otherwise 
independent of WIMP mass, the high mass DFs would 
be greater than the low mass DFs simply due to the 
prcfactor in Eq. (CI 3). In order to isolate the ef- 
fects of WIMP mass on the initial distribution of energy 
and angular momentum as well as subsequent rescatter- 
ing, we divide the throe DFs from the simulations with 
CTp^ = 10~^^ cm^ in Fig. 6(a) by and show these 
functions in Fig. 14. The low plateaus do not appear to 
be significantly different. There arc some discrepancies 
in the spikes, which are due to the low numbers of long- 
lived resonance-sticking WIMPs in each simulation. The 
high plateaus look relatively consistent with each other, 
given the large error bars. 



C. Maximum DF from Spin-Independent Solar 
Capture 

An important quantity to estimate is the maximum 
allowed DF consistent with experimental constraints on 
CTp ^. We expect the point in — a^^ yielding this max- 
imum DF to lie on the exclusion curve, but the maximal 
point is determined by the shape of the curve for the fol- 
lowing reason. The best limits on a^^ are shown in Fig. 4 
and come from the XENONIO (below = 40 GeV) and 
CDMS (above = 40 GeV) experiments [32, 64]. The 
exclusion curves reach a minimum of a^^ ~ 4x lO"*** cm'^ 
in the range m-^ = 20 — 70 GeV. Below these masses, the 
exclusion curve rises sharply due to kinematic reasons. 
Above m-^ 70 GeV, the exclusion curves rise oc 
because the flux of halo WIMPs at the detector goes as 

Since extensions to the standard model generically pre- 
dict > 100 GeV (with the notable exception of some 
gauge-mediated supersymmetry breaking models which 
predict ~ 1 keV [65-68]), we focus on this part of 
the exclusion curve [1, 2, 5]. In the previous sections, 
we found that (i) the high plateau dominates the DF 
at least up to a^^ = lO""*^ cm^, (ii) this plateau is a 
growing function of cross section until it reaches equi- 
librium for (jp^ > 10~*^ cm^, and (iii) for a fixed cross 
section, f(v)/n^ oc N^/n^, which reaches its maxinmm 
for ~ 2 TeV (Fig. 13(b)). If the CDMS exclusion 
curve in Fig. 4 were extended to higher mass, one would 
find that the exclusion curve hits = 10~^^ cm^ near 
= 2 TeV, which is exactly the point at which both 
equilibrium in the high plateau is achieved and N^/n^ 
reaches its maximum. 

In Fig. 15, we show the estimated DF for m-^ = 2 



TeV and = 10~^^ cm^, which we interpret as the 
maximum possible DF consistent with exclusion limits if 
spin-independent scattering dominates in the Sun. This 
DF is based on the DAMA simulation DF, appropriately 
scaled by WIMP cross section and mass. This DF yields 
a bound WIMP fraction (relative to the halo) which is a 
factor of ^ 4 greater than that of the DAMA simulation. 

In conclusion, we find that even for the maximal DF 
for bound WIMPs, the bound population is more than 
three orders of magnitude smaller than the total halo 
population at the Earth. 

D. Extension to Spin-Dependent Capture 

So far, we have only explored the dark matter DF 
in the case where WIMP-nucleon scatters in the Sun 
are dominated by spin-independent, scalar interactions. 
However, limits on the spin-dependent WIMP-proton 
cross section are 0(10'') times weaker than o-p ^, and spin- 
dependent cross sections arc generally higher than spin- 
independent cross sections in large parts of parameter 
space for well-motivated WIMPs. We showed that the 
low plateau of the solar capture DF, consisting of Jupiter- 
crossing WIMPs, grows as N^/n^ oc (7^^ or cr^^. Since 
the constraints on (7^^ are so much weaker than on cjp^, 
the low plateau could become large, reaching equilibrium 
when rescattering in the Sun occurs on timescales shorter 
than the time to pull the Jupiter-crossing WIMP perihe- 
lia outside of the Sun. 

In Fig. 16, we show a prediction for the low plateau 
for nix = 500 AMU and a^^^ = IQ-^^ cm^, if the only 
dependence of the DF on the WIMP cross section is 
f{v) oc Nq. The cross section is above the best m^ — ap^ 
unless > 1 TeV [69, 70], but is chosen to demonstrate 
an approximate maximum possible bound DF. At higher 
cross sections, the Sun becomes optically thick to WIMP 
scattering, at which point we expect the WIMP DF at 
the Earth to drop dramatically. The large central peak 
in the predicted DF arises from the nearly radial orbits. 
If the low plateau scales strictly with cross section un- 
til the Sun becomes optically thick, the Jupiter-crossing 
particles dominate the bound DF, and can swamp the 
unbound DF at low speeds {v < 50 km s~^). 

However, there are some indications within the simula- 
tions that the low plateau will grow less rapidly with cross 
section than in this simple model. Recall that ~ 98% 
of Jupiter-crossing WIMPs are ejected in the CDMS, 
Medium Mass, and Large Mass simulations, but a smaller 
fraction (w 73%) of WIMPs are ejected in the DAMA 
simulation. Therefore, a more careful estimate of the DF 
is required. 

To find how large the WIMP DF can get, we estimated 
the bound WIMP DF for various large spin-dependent 
cross sections {ffp^ > 10~^° cm^) using the DAMA sim- 
ulation as a starting point, since it has the highest 
and best statistics of all the spin-independent simula- 
tions. We scaled the total optical depth of each particle in 
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FIG. 13: In each plot, the red solid line denotes all species in the Sun, and the dotted blue line represents hydrogen, (a): The 
capture rate iV of WIMPs by the Sun for a'p^ = 10~*^ cm^, divided by the halo number density of WIMPs. The short solid 
black line gives the slope N /n^ oc m~^, the limiting slope for ^ niA for a nuclear species A. (6): The capture rate iV® to 
Earth-crossing orbits divided by the halo WIMP number density. 
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FIG. 14; DFs for the three simulations with ap' = 10""*^ cm^ 
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FIG. 15: The estimated maximum DF consistent with current 
exclusion limits if spin-independent scattering dominates in 
the Sun, for a WIMP mass = 2 TeV and WIMP-proton 



cross section ai. 
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the DAMA simulation by an estimate of the optical depth 
for a particular spin-dependent cross section. For parti- 
cles that were not on Jupiter-crossing orbits, we scaled 
the lifetimes by the ratio of the optical depth for the 
particular spin-dependent cross section and the DAMA 



optical depth. For the particles on Jupiter-crossing or- 
bits, we used the optical depth data from the DAMA 
simulation to find the approximate time at which each 
particle hit a total optical depth r = 1 for the new cross 
section, which we interpreted as the new WIMP lifetime. 
We calculated the DFs using the methods in Appendix 
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FIG. 16: Predicted geocentric DF if o-p 
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ing /(«) oc (jp^'"^^ for Jupiter-crossing orbits. This prediction 
is based on the output of the DAMA simulation. 



C, with the inclusion of a Monte Carlo treatment of the 
initial conditions to determine if captm'cd WIMPs scat- 
tered multiple times before they could leave the Sun. 

There are several assumptions in this approach. First, 
we used the initial distribution of semi-major axis and 
eccentricity derived from the DAMA simulation without 
any kinematic corrections due to the extreme mass dif- 
ference between hydrogen atoms and WIMPs. Thus, we 
tend to overestimate the Kozai contribution to the DF 
since scattering in the outer part of the Sun is suppressed 
for high . This also underestimates the contribution of 
Jupiter-crossing particles since the semi-major axis dis- 
tribution skews to higher a for large imbalances between 
the WIMP a. However, between = 60 AMU and 
= 500 AMU, the fraction of Earth-crossing particles 
that are also Jupiter-crossing only increases from 18.9% 
to 21.5% if the particles scatter only on hydrogen. 

Secondly, we did not recalculate optical depths for 
each passage through the Sun. This would be too time- 
consuming. Instead, we scaled the optical depths of each 
particle by the ratio of the scattering rate of i? = halo 
particles with the new cross section to the scattering rate 
of £^ = halo particles in the DAMA simulation. Since 
bound Earth-crossing particles do not have energies that 
vary significantly from E = relative to typical energies 
of unbound halo particles, using the ratio of the scat- 
tering rates to scale the DAMA optical depths should 
be a reasonable proxy for finding optical depths for spe- 
cific paths through the Sun. However, this approximation 
does neglect any differences in the radial distributions of 
hydrogen and heavier elements in the Sun, as well as any 
kinematic effects due to scattering off hydrogen rather 
than heavier atoms. 
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since the optical depth of the Sun for (jp = 10 
approximately the same as cTp ^ = 1.3 x 10""^^ cm^ 



We estimated DFs for = 60 AMU at a^^ = 
1.3 X 10-39,10-38,10-37, and 10'^^ cm^, and then ex- 
trapolate the results to other WIMP masses by rescaling 
the DFs by N§ ("ix)' ^^^^ scattering halo WIMPs 
on hydrogen to reach bound. Earth-crossing orbits. The 
cross section (7^^ = 1.3 x IQ-^^ cm^ yields similar same 
optical depths in the Sun as a^^ = IQ-*^ cm^. We used 
50 bootstrap resamplings for each spin-dependent cross 
section to estimate the DFs. The results are shown in 
Fig. 17, displaying f{v)/n^ for each cross section against 
the geocentric unbound DF. 

There are several key points this figure. The central 
part of the DF for each cross section (t; = 30 — 45 km 
s-^) is approximately independent of cross section, which 
is what one would expect if Kozai cycles dominate this 
region and particles have lifetimes of at least one Kozai 
cycle. This region is relatively unaffected by multiple 
scatters before the WIMPs exit the Sun for the first time 
because the particles on Kozai cycles originate in a part 
of the Sun that still has very low optical depth, even for 
the highest cross section considered. The peak near 50 
km s-^ is due to nearly radial Jupiter-crossing orbits. 
The spikes in the low plateau grow for a while and then 
disappear, a consequence of rescattering in the Sun before 
WIMPs can stick to resonances. 

The most striking result of Fig. 17 is that the low 
plateau is quite a bit lower than the naive prediction 
in Fig. 16. It appears that, while the low plateau 
does rise for large WIMP-proton cross sections, rescat- 
tering in the Sun plays an integral role in severely reduc- 
ing Jupiter-crossing particle lifetimes. We find that the 
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low plateau reaches approximately its maximum value 
if cTp^ = 10~^^ cm^. Even though the low plateau is 
still very slowly increasing between (7^ ^ = lO"^'^ cm^ for 

V < 50 km s~^, the plateau actually decreases between 
going from (7^^ = lO^'^"^ cm^ to a^^ = 10"^^ cm^. This 
is because WIMPs with geocentric speeds v > 50 cm 
are retrograde with respect to the planets in the solar sys- 
tem, and the torques from the planets are less effective 
for retrograde than prograde WIMPs. Thus, the time for 
WIMP perihelia to exit the Sun is longer for retrograde 
WIMPs than prograde WIMPs, and so the probability 
for a retrograde WIMP to rescatter in the Sun before its 
perihelion exits the Sun for the first time is significantly 
higher than for a prograde WIMP. Therefore, the maxi- 
mum low plateau occurs for a^^ ~ 10~^^ cm^, or about 
aSi ^ 10-38 cm2. 

Combining these results with the maximum DF for 
spin- independent solar capture in Fig. 15, we find that 
particles captured to the solar system by elastic scatter- 
ing in the Sun are only small population relative to the 
halo population at the Earth, even if the spin-dependent 
WIMP-proton elastic scattering cross section is quite 
large. Improving on the approximations we used in this 
section is unlikely to change this conclusion. 



V. THE DIRECT DETECTION SIGNAL 

Direct detection experiments look for nuclear recoil of 
rare WIMP-nuclear interactions in the experimental tar- 
get mass. The WIMP-nucleus scattering rate per kg of 
detector mass per unit recoil energy Q can be expressed 
as [cf. 1] 

where da^/dQ is the differential interaction cross section 
between a WIMP and a nucleus of mass niA and atomic 
number A, and v is the velocity of the dark matter par- 
ticle with respect to the experiment. The lower limit to 
the integral in Eq. (78) is set to 

Vrmn = (m AQ /2n\y/^ , (79) 

the minimum WIMP speed that can yield a nuclear recoil 
Q, The dark matter DF at the Earth is /(x, v). 

In this section, we will determine the maximum possi- 
ble contribution of the bound DF to the direct detection 
rate. We focus on the maximum event rate from bound 
WIMPs instead of exploring how the bound WIMP event 
rate depends on WIMP mass and elastic scattering cross 
section since we expect the event rate to be small. We are 
interested in both the total excess signal due to bound 
WIMPs for particular experiments, as well as the con- 
tribution to the differential event rate, since the latter is 
important for determining the WIMP mass [31, 71]. 



We focus on directionally-insensitive direct detection 
rates for spin-independent interactions, but the results 
of this section can be applied qualitatively to spin- 
dependent interactions as well. There is another class 
of direct detection experiment that is directionally sen- 
sitive [72-76]. In principle, the bound WIMPs should 
leave a unique signal in such experiments (see Fig. 7), 
but it would be challenging to measure this given the 
small bound WIMP density, current errors in directional 
reconstruction, and high energy thresholds. 

We calculate the bound WIMP event rate for = 
500 AMU and a^^ = 10~^^ cm^ and a high spin- 
dependent proton cross section {(7^^ = lO^''^ cm^, ap- 
proximately the point at which the Sun becomes optically 
thick to WIMPs). We choose this point in parameter 
space because it yields the largest DF due to WIMPs 
bound by solar capture. The event rate can simply be 
scaled for lower (or higher) spin-independent cross sec- 
tions. The scaling for other values of and (7^^ is 
different, but can be easily determined. 

The geocentric bound DF is anisotropic. Therefore, 
to translate the DF outside the sphere of influence of 
the Earth to the corresponding DF at the detector, one 
should use the mapping technique in Appendix C, aver- 
aged over the detector's daily motion about the Earth's 
rotation axis. However, using the isotropic mapping in- 
stead of the full six-dimensional mapping in Appendix 
C produces errors in dR/dQ of at most a few percent. 
Therefore, we use this simplification for the bound WIMP 
DF at the surface of the Earth in calculating dR/dQ. 

In Fig. 18, we show the maximal direct detection signal 
due to solar captured WIMPs if = 500 AMU (lower 
two lines). We find direct detection rates assuming ^^^Xe 
and ^^Ge targets, since the current and planned experi- 
ments most sensitive to the spin-independent (and spin- 
dependent neutron) cross section have isotopes of either 
Xe or Ge as their target mass. For comparison, we also 
plot the event rate expected for the halo DF in Eq. (49). 
We find that bound WIMPs can only enhance the direct 
detection rate at very small Q, and that the enhancement 
is largest at the smallest recoil energies. For both the ger- 
manium and xenon targets, the maximum enhancement 
to the total event rate is ~ 0.5% at Q = 0. This enhance- 
ment is actually disproportionally large compared to the 
enhancement in the local WIMP number density due to 
bound WIMPs, which is n^^bound ~ ^O~'^n^jiaioi since 
incoherence in the WIMP-nucleon interaction for large 
nuclei suppresses the elastic scattering cross section for 
high speed halo WIMPs. We also show the experimental 
analysis windows for the recent XENONIO and CDMS 
experiments in this figure [32, 64]. The current analysis 
threshold of the CDMS experiment is too high to de- 
tect bound WIMPs. If this experiment and its successor 
SuperCDMS could push down their analysis thresholds, 
as other germanium-based rare event experiments have 
(e.g., CoGeNT [33]), bound WIMPs may be observed. 
At Q = 4.5 keV, the current analysis threshold for the 
XENONIO experiment, the boost to the differential event 
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FIG. 18: The differential direct detection rate for = 500 
10"*'^ cm^ assuming the DF is dominated by 



AMU and a^' 

spin-dependent scattering in the Sun with cjp^ = 10-'^^ cm^. 
The shaded region indicates the XENONIO analysis region 
[32], and the vertical dashed line indicates the lower limit to 
the CDMS analysis window (which extends to Q = 100 keV) 
[64]. 



rate is ^ 0.1%, and the total boost in their analysis win- 
dow is ~ 10""^%. Thus, the bound particles only neg- 
ligibly increase the total event rate (integrating dR/dQ 
over the range of Q's allowed in the analysis window), if 
at all. Estimates of the WIMP mass and cross section 
from direct detection experiments will not be affected by 
solar captured particles. 



VI. THE NEUTRINO SIGNAL FROM WIMP 
ANNIHILATION IN THE EARTH 

WlMPs may accumulate and annihilate in the Earth. 
The signature of WIMP annihilation will be GeV to TeV 
neutrinos. Neutrino observatories (e.g., Antares [77], Ice- 
Cube [78]) are sensitive to the Cerenkov radiation of 
muons created in charged-current interactions of muon 
neutrinos in and around the experiment. 

The neutrino flux at a detector on the surface of 
the Earth is proportional to the annihilation rate F of 
WIMPs trapped in the Earth. Finding F requires solv- 
ing a differential equation for the number of WIMPs N 
in the Earth, described by 

N ^C-2T, (80) 

where the capture rate of WIMPs in the Earth by elastic 



C= / d^x 



l'/<l'csc(x) 



dVdf^^^n^(x)z 



X /(x,v,i). (81) 



Here, da^/dri is the WlMP-nucleus elastic scattering 
cross section for nuclear species A and v is the relative 
speed between the WIMP and a nucleus. The number 
density of species A is described by riyi(x). The cutoff 
in the velocity integral reflects the fact that the WIMP's 
speed after scattering Vf must be less than the local es- 
cape velocity Wesc(x). 

If the WIMP DF is time-independent, the annihilation 
rate goes as 



F= ictanh2(tAe), 



where 



te 



(CCa) 



-1/2 



(82) 



(83) 



is the equilibrium timescale and Cq is a constant that 
depends on the distribution of WIMPs in the Earth and 
is proportional to the annihilation cross section. 

While the contribution of bound particles to the direct 
detection rate is expected to be minuscule, it is not unrea- 
sonable to expect that the bound particles could notice- 
ably boost the neutrino-induced muon event rate from 
WIMP annihilation in the Earth. Because the Earth's 
gravitational potential is shallow, it is difficult for halo 
WIMPs to lose enough energy during collisions with the 
Earth's nuclei to become bound unless the WIMP mass 
is nearly equal to the mass of one of the abundant nu- 
clear species in the Earth [79]. For WIMPs with mass 
TO^ > 400 GeV, only WIMPs bound to the solar system 
may be captured in the Earth. 

In Fig. 19, we show the capture rate (Eq. 81) of 
WIMPs in the Earth as a function of mass for = 
10"**^ cm^ for several different WIMP DFs. We use the 
potential and isotope distributions in Encyclopaedia Bri- 
tannica [80] and McDonough [81]. The lowest line shows 
the capture rate of only the unbound WIMPs in the so- 
lar system. The peaks in the capture rate correspond 
to points at which the WIMP mass is nearly exactly the 
same as a one of the common elements in the Earth, 
of which the iron peak {mpe « 56 AMU = 53 GeV) 
is especially prominent. The long-dashed line represents 
the capture rate for both unbound particles and particles 
bound to the solar system by spin-independent scatter- 
ing in the Sun. We show extrapolations to the regime 
in which spin-dependent scattering dominates in the Sun 
with the short dash-dotted and long dash-dotted lines, 
representing a^^ = 1.3 x 10^^^ cm^ and cr^^ = 10^^^ 
cm^ respectively. We included unbound WIMPs in those 
estimates. 

From Fig. 17, we note that the high plateaus in the 
DF's if (Tp^ > 10^"^^ cm^ are nearly identical; the main 
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FIG. 19: Capture rate of WIMPs in the Earth as a function of 
WIMP mass for ap' = 10~*^ cm^. All capture rates include 
the capture of unbound halo WIMPs as well as the capture 
of bound WIMPs. 



reason for the difference in the capture rate is the low 
speed DF of Jupiter-crossing WIMPs. In fact, the cap- 
ture rate is extremely sensitive to the DF of the lowest 
speed WIMPs. For the relatively low capture rates in 
Fig. 19, ie > t©, so r cx C^. Even small variations in the 
low speed WIMP DF can lead to large variations in the 
event rate at a neutrino telescope. 

To estimate a plausible range of muon event rates given 
the capture rates in Fig. 19, we explore part of the MSSM 
parameter space. We can in principle explore other mod- 
els, but the MSSM yields, on average, somewhat larger 
spin-independent cross sections. Given that iron is the 
most common element in the core of the Earth, and oxy- 
gen, silicon, and magnesium the most common element 
in the mantle, none of which has spin-dependent inter- 
actions with WIMPs, only in WIMP models with appre- 
ciable spin-independent interactions will capture in the 
Earth be relevant. 

We scan MSSM parameter space to estimate the 
neutrino- induced muon event rate for neutrino telescopes 
from neutralino annihilation in the Earth using routines 
from the publicly available DarkSUSY v.5.0.2 code [82]. 
The code can also check whether a model described by a 
set of SUSY parameters is consistent with current collider 
constraints and relic density measurements. We describe 
our scans in more detail in Paper III. 

To estimate the muon event rate in a neutrino tele- 
scope, we set the muon energy threshold to EjJ^ — 1 
GeV. This is somewhat optimistic for the IceCube exper- 
iment [36, 83] unless muon trajectories lie near and ex- 



actly parallel to the PMT strings, but it is reasonable for 
the more densely packed water experiments (e.g., Super- 
Kamiokande). The signal drops sharply with increasing 
muon energy threshold [84]. We assume that the mate- 
rial both in and surrounding the detector volume is either 
water or ice, since the largest current and upcoming neu- 
trino telescopes are immersed in oceans or the Antarctic 
ice cap. We include all muons oriented within a 30° cone 
relative to the direction of the center of the Earth. 

In the following figures, we present muon event rates 
in neutrino telescopes for various DFs. In Fig. 20, we 
show the event rates for WIMPs unbound to the solar 
system. The solid black line on Fig. 20 represents the 
most optimistic flux threshold for IceGube [36, and ref- 
erences therein] . To show how the event rates depend on 
the SUSY models for a given spin-independent cross sec- 
tion, we mark the models on the figure according to which 
direct detection experiments bracket the cross section for 
a given neutralino mass. The open circles correspond to 
SUSY models with above the that lie above the 2006 
GDMS limit [85], which is shown in Fig. 4. The triangles 
are models for which a^^ lies between the 2006 GDMS 
limit and the current best limits on tip ^ (a combination of 
XENONIO [32] and GDMS [64] hmits), and squares de- 
note models consistent with all current direct detection 
experiments. 

It appears that no halo WIMPs from any of the mod- 
els found in our scan of the MSSM consistent with exper- 
iments would produce an identifiable signal in IceGube. 
We cannot say that it is impossible for neutralino WIMPs 
from the halo to be observed by IceGube or other km"^- 
scale experiments, since we are only sampling a small part 
of the vast SUSY parameter space, but Fig. 20 suggests 
that it is not likely. 

In Fig. 21, we show the muon flux for WIMPs cap- 
tured in the Earth from the halo or from the population 
of bound WIMPs. We calculate the muon event rate 
with the bound DF for = 1.3 x 10"^^ cm^ no mat- 
ter what the actual tXp in the model is since this is near 
the maximum spin-dependent cross section found in the 
parameter scans, a^^ is almost always small enough that 
the DF due to spin-independent scattering in the Sun is 
subdominant to the spin-dependent-derived DF. There- 
fore, the points in Fig. 21 are almost entirely upper lim- 
its to the solar captured WIMP event rates. This figure 
is almost indistinguishable from Fig. 20. We find that 
the maximum enhancement over the halo WIMP event 
rate is of order 20%. Thus, the solar captured WIMPs 
produce almost no enhancement in the neutrino-induced 
muon event rate. 

One caveat to this pessimistic result is that we esti- 
mated the event rate using only the fiux of muons from 
outside the detector volume. However, Bergstrom et al. 
[84] suggest that muons created inside the detector vol- 
ume may dominate the signal for smaller WIMP masses 
{m^ < 300 GeV) in large (km'^) telescopes, although the 
exact enhancement has not been calculated. But, the en- 
hancement of the event rate due to bound WIMPs over 



28 



CO 



1000 



100 



10 



0.1 



0.01 



CO 
K 

> 



o 

^ 0.001 



0.0001 



Unbound 



" "\ 

■■ ■ ■■ l^' ■ ■ 

I # . ■ ■ 

t^:...ic:..-.i 



100 



Mass [GeV] 



1000 



FIG. 20: Muon event rates from halo WIMPs unbound to the 
solar system. Open circles mark MSSM models for which ap' 
is above the 2006 CDMS limit [85], filled triangles mark those 
with limits between that limit and the current best limits on 
a^' (set by XENON 10 for < 40 GeV [32] and CDMS for 
> 40 GeV [64]), and filled squares denote models consis- 
tent with the best fimits on elastic scattering cross sections. 
The solid line is an optimistic detection threshold for the Ice- 
Cube experiment [36, and references therein]. 
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halo WIMPs will be fixed and small. 



VII. DISCUSSION 



A. Comparison with Damour and Krauss 

Here, we compare the simulation results with the semi- 
analytic predictions in Damour and Krauss [37]. 

Damour & Krauss neglected the population of WIMPs 
on Jupiter-crossing orbits, arguing that it would be short- 
lived because of the strong perturbations from Jupiter. 
This argument is plausible, but it neglects the impor- 
tance of long-lived WIMPs on resonant orbits. The pres- 
ence of long-lived WIMPs on resonances suggests that 
Jupiter-crossing WIMPs may be important for a^^ ~ 
IQ-^i cm2 (o-^^ - 10^39 cm^). However, such WIMPs 
are unlikely to contribute significantly for much larger 
or much smaller cross sections. For much larger cross 
sections, long-lived WIMPs should be exceedingly rare; 
they are likely to rescatter and thermalize in the Sun be- 
fore Jupiter can pull the perihelia out of the Sun. For 
smaller cross sections, the rate of scattering of WIMPs 
onto Jupiter-crossing orbits is negligible. 

Before we describe where our results diverge from 
Damour & Krauss for a < a%/2 k, 2.6 AU, we reem- 
phasize the main points of their work. They found that 
the main enhancement to bound WIMP DF came from 
a small fraction, ~ 0.1%, of WIMPs scattered onto or- 
bits with 0.5 AU < a < 2.6 AU on Kozai cycles. They 
assumed that these WIMPs, which originated in the out- 
skirts of the Sun, had lifetimes at least as long as the 
age of the solar system ~ 4.5 Gyr. For this range of 
semi-major axes, we found two major differences between 
their work and ours. 

First, we find that WIMPs with 1.5 AU < a < 2.6 AU 
are not well described by pure Kozai cycles due to sig- 
nificant interactions with mean-motion resonances. Un- 
less the WIMP-nucleon cross section is large {(7^^ ^ 
10^^^ cm^, allowed if < 5 GeV or 



> 10 TeV 

[32, 64]; a^" » IQ-^^ cm^), most of the WIMPs in this 
semi-major axis range have lifetimes ~ 100 times longer 
than if the Sun were an isolated body. However, this still 
does not increase the DF at the Earth as much as if the 
1% of WIMPs in this semi-major axis band (the fraction 
of WIMPs initially scattered onto 1.5 AU < a < 2.6 AU 
which were on Kozai cycles in [37]; a higher fraction of 
large semi-major axis WIMPs are on Kozai cycles than 
WIMPs with lower semi-major axis) had lifetimes ex- 
tending to Iq. 

To show why, we use the following argument. The 
increase in the number density of WIMPs at the Earth 
n' over the number density without tidal torques n is 
roughly described by 



FIG. 21: Muon event rates including bound WIMPs. Symbols 
mark the same models as in Fig. 20. 



n 
n 



(84) 



where e/ is the fraction of WIMPs disturbed enough 
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from their orbits to have significantly longer lifetimes 
in the solar system. The factor Et describes the 
increase in the WIMP lifetime. Typically, Et w 
i0)/min(imed,i0), where Uned is the median 
lifetime of the WIMPs in the absence of gravitational 
torques and t'^gd is the median lifetime with gravita- 
tional torques. For our simulations, ^ 1 since most 
WIMPs with 1.5 AU < a < 2.6 AU were on quasi-Kozai 
orbits and Et ~ 100, implying that n'/n ^ 100. How- 
ever, the Damour and Krauss [37] prediction would be 
e/ ~ 0.01, and Et - (to/llO^ yr) = 4.5 x 10^, implying 
that n'/n ~ 4.5 x 10'', a factor of ^ 500 larger than what 
we found in our simulations. 

However, the main reason that the density of bound 
WIMPs is much smaller than estimated by Damour and 
Krauss is that particles with a < 1.5 AU on Kozai cy- 
cles have lifetimes that are much less than the age of the 
solar system. This is due to the fact that the typical in- 
tegrated optical depth per Kozai cycle is non-negligible, 
so a WIMP undergoes only a finite number of Kozai cy- 
cles before rescattering in the Sun. There are two im- 
portant timescales relevant to estimating the lifetimes of 
WIMPs on Kozai cycles for a given WIMP-nucleon scat- 
tering cross section. 

First, in the point mass three-body problem, the period 
of Kozai cycles are of order [cf. 86] 



Toe 



P M' 



(85) 



Here, P denotes the orbital period of a particle and Ptj, 
represents the orbital period of Jupiter. For typical par- 
ticles, T < 10^ yr. 

The other important timcscalc is the timescale on 
which the orbital perihelion is moved out of the Sun. 
Although the optical depth in the outskirts of the Sun is 
extremely low (r ~ 10~^ for an orbit with w O.7i?0, 
r ^ 10~^ for Tp K. Q.QRq in the DAMA simulation, and 
even lower in the other simulations), it takes many orbital 
periods for Jupiter to pull the perihelia out of the Sun, 
hence making the optical depth per Kozai cycle much 
larger than the optical depth for a single passage through 
the Sun. 

The rate of change in the angular momentum of a 
WIMP is 



(86) 



where K\ is the torque on the particle orbit by Jupiter. 
The torque is larger at aphelion Ta for particles with a < 
a\/2 than at any other point in the orbit, so the average 
torque can be approximated by its value at aphelion 



(87) 
(88) 



applied at aphelion, where we have expanded the poten- 
tial to the / = 2 term in the spherical harmonic expansion 



(Eq. 41). The angular momentum must change by of or- 
der 



(89) 



for perihelia to be external to the Sim. In reality, since 
WIMPs on Kozai cycles originate at distances from the 
center of the Sun r > O.5R0, Eq. (89) should have a 
small (^ 0.1 — 1) coefficient in front. Therefore, if the 
torques are coherent, the total time it takes for a WIMP 
to have its first perihelion outside the Sun is 



At . 



AJ 



(90) 



Using the expressions for K\ and AJ in Eqs. (88) and 
(89), we find 



Ai Mq f a \ "^/^ / a 
~ 10'', for a = 1 AU. 



(91) 
(92) 



Thus, a particle passes through the Sun many times dur- 
ing each Kozai cycle. In the simulations, we find that 
the total optical depth per Kozai cycle is ^ 10^ — 10'' 
times the optical depth at maximum eccentricity. Even 
if the optical depth at maximum eccentricity is only 
10"*" - 10"^ per orbital period (typical of the DAMA 
simulation), the total optical depth per Kozai cycle is 
^ 10"'^. It only takes about 1000 Kozai cycles for such 
a particle to rescatter in the Sun. The result is that the 
lifetimes of particles are typically less than the age of the 
solar system (~ 100 Myr), and as such cross the Earth's 
orbit a factor of ~ 50 times than predicted by Damour 
and Krauss. 

To compare our results to Damour & Krauss, we use 
Eq. (84). They find that e/ ~ lO'^ of WIMPs with 
0.5 AU < a < 1.5 AU initially captured in the Sun will 
be on a Kozai cycle. For their typical WIMP-proton cross 



section a^^ 



10- 



so tmed ~ lO"* yr. 



and Et ~ 4.5 x 10^. Thus, Damour and Krauss expect 

n'/n - 10^ - 10-*. 



However, for the same cross section, we find <' 



10« 



yr, such that Ei ^ 10^. Thus, n'/n ^ 100, which is ap- 
proximately the upper limit of what is found in the simu- 
lations. In general, we find n'/n somewhat smaller than 
~ 100, both because Et decreases as ap moves farther 
below the equilibrium value {af,^ ~ 10"''^ cm^), and be- 
cause the median lifetime of WIMPs not on Kozai cycles 
but drawn from the same a and rp as the Kozai cycle 
WIMPs is a bit higher than the population of WIMPs 
with a < oi,. as a whole. 

We find that we can recover the Damour and Krauss 
[37] estimates of the maximum increase in direct detec- 
tion experiments if the Kozai WIMPs in our simulations 
had never scattered. For the DAMA simulation, the 
median Kozai WIMP lifetime is just short of 100 Myr 
(Fig. 10). If these WIMPs had instead rescattered on 
timescales longer than the age of the solar system, then 
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wc would expect the DF to have been larger by a factor 
of ~ 50 — 100. We found that the maximum increase to 
the differential direct detection rate dR/dQ (Eq. 78) was 
~ 0.5% of the halo event rate. If the DF were larger by 
this factor of 50-100, then the bound WIMPs would add 
an additional 25 — 50% of the halo event rate at small Q, 
consistent with what is found by Damour & Krauss. 

We can also recover the large neutrino event rate from 
WIMP annihilation in the Earth found by Bergstrom 
et al. [38] using the Damour & Krauss results. We 
found that for MSSM models consistent with limits on 
the WIMP-nucleon elastic scattering cross section, the 
capture rate of solar-captured WIMPs in the Earth was 
a maximum of about 10% that of the halo, for w 100 
GcV. If the DF were a factor of 50— 100 higher, the solar- 
bound WIMP capture rate would be 5 — 10 times higher 
than the halo capture rate. Since tg > (Eq. 83) for 
such capture rates, the annihilation rate of WIMPs in the 
Earth would scale as F oc C^, leading to an increase in 
the neutrino flux in neutrino detectors of 25 — 100 times 
the halo event rate, consistent with what was found by 
Bergstrom et al. However, we note that even if the en- 
hancement were that high. Fig. 21 shows that this signal 
would fall below the IceCube flux threshold for WIMP 
models consistent with experimental constraints. 



solar system, which we show in more detail in Chapter 
5 of [52]. However, encounters with other planets can 
push geocentric WIMP speeds below t; = 10 km s~^ by 
increasing J^. While the presence of a tail in the DF at 
low geocentric speeds is not significant for direct detec- 
tion event rates, it can have a disproportionate eff'ect on 
the capture rate of WIMPs in the Earth. 

Secondly, the overall number density of bound WIMPs 
may change, depending largely on how efficient the plan- 
ets are at increasing (or decreasing) the lifetimes of 
WIMPs in the solar system (Eq. 84). We will argue 
below that the true number density of WIMPs at the 
Earth is unlikely to be much larger or smaller (within 
factors of a few) than that estimated from simulations 
using a toy solar system. 

We divide the discussion into three parts: (i) WIMPs 
with initial < 1.5 AU, (ii) WIMPs with 1.5 AU < 
tti < 2.6 AU (quasi-Kozai WIMPs in the toy solar sys- 
tem), and (iii) Jupiter-crossing WIMPs. Without further 
simulations, though, it is not possible to tell exactly by 
how much the DF will change. Hence, we also discuss 
the challenges involved in simulating WIMPs in a more 
realistic solar system. 

1. a < 1.5 AU 



B. Planets 

Of course, all of the conclusions in this work are based 
on simulations in a toy solar system, consisting of Jupiter 
on a circular orbit about the Sun. Dynamics in the solar 
system are much more complex, both because Jupiter has 
non-zero eccentricity and inclination and because other 
planets arc present. Bodies may have close encounters 
with any planet within its aphelion, and may be influ- 
enced by additional mean-motion and secular resonances 
[e.g., 62, 87-90]. The combination of these effects yields 
far greater diversity of orbits in the real solar system than 
what we found in the toy solar system. 

There are two qualitatively different ways in which a 
more realistic treatment of the solar system could change 
the WIMP distribution at the Earth. First, additional 
parts of phase space become accessible. While it is a tri- 
umph of our numerical methods that the Jacobi constant 
is conserved to high accuracy in our simulations, the con- 
servation of the Jacobi constant restricts the range of mo- 
tion for WIMPs. For example, WIMPs with a < a%/2 ex- 
perienced only minor fluctuations in the semi-major axis 
because they never encountered Jupiter closely enough to 
experience large energy changes. Thus, according to the 
definition of the Jacobi constant, Eq. 20, even WIMPs on 
quasi-Kozai cycles only experienced relatively minor per- 
turbations to Jz- This meant that WIMPs not crossing 
Jupiter's orbit had heliocentric velocities perpendicular 
to the Earth's motion, restricting the geocentric speeds 
V ^ 30 km . Jupiter-crossing WIMPs were re- 

stricted to geocentric speeds t; > 10 km s~^ in the toy 



The DF of solar-captured WIMPs could be greatly in- 
creased if the planets other than Jupiter were to either 
(i) pull a larger percentage of particles out of the rescat- 
tering peak and onto orbits that only occasionally enter 
the Sun or (ii) extend the lifetimes of particles that al- 
ready did exit the Sun in the toy solar system simulations. 
Here, we discuss three mechanisms for pulling additional 
WIMPs out of the Sun: (i) close encounters with inner 
planets, (ii) changes to the Kozai structure by other plan- 
ets, and (iii) additional secular resonances. Then, we will 
estimate the lifetimes of such WIMPs. 

Close encounters: Here, we describe how random-walk 
encounters with planets can pull WIMPs that were in the 
rescattering peak in our simulations onto long lifetime 
orbits in the solar system. Close encounters with the 
inner planets can alter the WIMP angular momentum 
with respect to the Sun. Ignoring resonant phenomena 
in the solar system, the close encounters can be treated 
as a diffusion problem. We use the rms change in an- 
gular momentum as a function of time to estimate the 
timescales on which WIMP perihelia are pulled out of the 
Sun. Modeling WIMP-planet interactions as two-body 
encounters, each time a WIMP of heliocentric speed v 
crosses a planet's orbit, the WIMP's planet-centric speed 
u changes in the direction perpendicular to u by 



where b is the impact parameter. Since WIMPs with 
a < 1.5 AU are on ex tremely eccentric orbits, to good 
approximation, u = ^/v"^ + Vp, where Vp = GMq/op. 
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The change in planet-centric speed can be related to the 
change in heliocentric speed by 
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As a rough approximation, the change in angular mo- 
mentum per encounter is thus 



6J ~ ap6v. 



(96) 



We use the approximation that the angular momentum 
undergoes a random walk to estimate the timcscale on 
which a particle's angular momentum changes by of order 
AJq ~ ^JGMqRq (Eq. 89) in order for the orbital 
perihelion to lie outside the Sun. The rms change in 
angular momentum will go as 



((AJ)2) ~ N{5Jf 



(97) 



where the particle encounters planet P with an impact 
parameter b or less a total of N times in a time span t. 
In general, 
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where P-^^ is the orbital period of the WIMP. The factor 
{b/ap)^ is the probability per WIMP period that the 
WIMP comes within a distance b of the planet. Thus, 
with some rearranging, we find 
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where the factor of 10 comes from the heretofore ig- 
nored Coulomb logarithm (see [91]). The singularity at 
a = ap/2 is artificial and would vanish in a more care- 
ful treatment of WIMP-planet encounters. Thus, the 
timescale for WIMPs to diffuse out of the Sun due to 
the action of planet P is 



id/yr ~ 



\Mp) \ap ) Vae 
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x(2-^). (100) 



For both the Earth and Venus, Mp/Mq ~ 3 x 10^ 
and R^/ap ^ 10~^, yielding a diffusion time td ^ 10^ 
yr for a ~ 1 AU. The timcscalcs for Mercury and Mars 
are ta ~ 10"^^ and ~ 10^° years respectively. Thus, an- 
gular momentum diffusion is dominated by the Earth 
and Venus. To estimate the impact on the number 
density, we must find e/, the fraction of WIMPs with 
a < 1.5 AU that may be perturbed out of the Sun. For 



10^ yr. 
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implying that 

2 ....10-3. 



the DAMA simulation 

e/ ~ tmed/td ~ 10"^. For a^^ = 10"'" cm', e/ 

To estimate the impact of this population on the num- 
ber density of bound WIMPs, wc must also estimate Et, 
the ratio of the median lifetime including the gravita- 
tional effects of the planets to the lifetime if the Sun 
were isolated. Still ignoring resonances, we estimate the 
rms timescale for WIMPs to be ejected from the solar 
system once the perihelia are outside the Sun, 
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Since 
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we use the expression for Sv in Eq. (95) to find 
Using the expression for N in Eq. (98), we find that 
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where again wc have included a factor of 10 for the 
Coulomb logarithm. The inner planets which will per- 
turb the orbits the most are Venus and the Earth, yield- 
ing ejection timcscales of tej ~ 10"'^'^ yr, longer than the 
age of the solar system. This yields Et ^ a few x 10^ 
for the DAMA simulation and Et ^ a few x 10** if 
(Tp^ = 10~^^ cm^. Combined, this would yield n'/n <^ 
a few x 10, where n is the number density of the rescat- 
tering peak WIMPs in the toy solar system simulations. 
This is of the same order as the increase in the bound 
WIMP DF due to Kozai cycles in our simulations. 

However, there arc reasons to believe that Et is in 
fact significantly smaller than these estimates suggest. 
First, if a WIMP can diffuse out of the Sun, it can 
also diffuse back in. Secondly, once a WIMP becomes 
Jupiter-crossing, it will be ejected from the solar system 
on timescales of ~ Myr, which is essentially instanta- 
neous. 

Thirdly, studies of Near Earth Object (NEO) orbits 
show that once small bodies reach o > 2 AU, they 
are driven into the Sun on rather short timescales, 
1 — 10 Myr, mostly by secular resonances but also by 
mean- motion and Kozai resonances [87, 89]. Given that 
WIMPs have significantly higher eccentricity that the 
typical NEO, the timescale to drive a WIMP back into 
the Sun via resonances may be shorter. On the other 
hand, such WIMP orbits have high speeds relative to the 
planets, while the low eccentricity, prograde, low inclina- 
tion NEO orbits have relatively low speeds. Hence, NEOs 
will be more efficiently gravitationally scattered onto the 
mean-motion and secular resonances that drive up the 
eccentricity. In spite of this latter effect, it is likely that 
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WIMPs will be scattered back into the Sun on timescales 
shorter than the age of the solar system. If the integrated 
optical depth in each instance that the WIMP perihelion 
is driven into the Sun (i.e., that the WIMP experiences 
many Sun-penetrating orbits each time a resonance ini- 
tially drives the WIMP into the Sun), the lifetime of the 
WIMPs will be less than the age of the solar system, 
hence reducing Et- 

Fourthly, Gladman et al. [89] have identified additional 
resonances that drive some NEOs of a < 1.9 AU into the 
Sun without first boosting the semi-major axis above a = 
2 AU. This will reduce Et for WIMPs with a < 1.9 AU. 

However, WIMPs can survive many passes through the 
Sun before scattering with solar nuclei onto uninterest- 
ing orbits. The timescale for rescattering in the Sun de- 
pends crucially on how many passages WIMPs can make 
through the Sun before gravitational torques from the 
planets pull the perihelia out again. 

In general, it appears that the lifetimes of WIMPs with 
o < 1.5 AU initially pulled out of the Sun by angular mo- 
mentum diffusion will be shorter than those predicted by 
arguments based on energy diffusion, although quanti- 
fying this is difficult without a full solar system Monte 
Carlo simulation. Even if WIMP lifetimes were domi- 
nated by diffusion instead of the effects listed above, the 
boost to the DF would only just be comparable to that 
due to Kozai cycles in the toy solar system. 

Changes to the Kozai structure: Next, we consider 
changes to the Kozai structure caused by planets other 
than Jupiter. Both the inner and outer planets can affect 
the structure of Kozai cycles. However, torques from the 
outer planets other than Jupiter are unlikely to change 
the number of particles whose perihelia exit the Sun. As 
demonstrated in Eq. (88), the torque on a particle by a 
faraway planet goes as oc Mpo^a^^ , where Mp and 
ap are the mass and semi-major axis of the planet, and 
a is the semi-major axis of the particle. A planet will 
provide a torque 

relative to the torque from Jupiter. Even Saturn, the 
next largest planet in the solar system, and the second 
nearest gas giant to the Earth, will only produce a torque 
about 5% that from Jupiter. Jupiter dominates the tidal 
field for particles that do not cross the orbits of the outer 
planets, and so it dominates the structure of the Kozai 
cycles. 

Among the inner planets, Michel and Thomas [44] find 
that the Earth and Venus can dominate the structure of 
the Kozai cycles if the semi-major axis of the particle 
is near the semi-major axis of either planet, the initial 
eccentricity of the particle orbits is small, and the max- 
imum inclination of the orbit is low. However, WIMPs 
tend to have power-law distributed semi-major axes, high 
eccentricities, and are scattered isotropically in the Sun. 
Therefore, we expect that the extra planets will not in- 



crease the number of particles on Kozai cycles in the inner 
solar system. 

Secular resonances: There are additional secular reso- 
nances in the full solar system that do not appear in the 
circular restricted three-body problem considered in this 
work. These occur when the rate of change of either the 

longitude of perihelion {zu) or of the longitude of the as- 
cending node {(l) of the WIMP is almost equal to that of 
one of the planets. The evolution of NEOs is greatly af- 
fected by the secular resonances with Jiipitcr and Saturn, 
although several authors show that other resonances are 
also important [87, 89, 92-94]. There are complications 
in interpreting and extending results from NEO simula- 
tions. For example, most analytic and numerical effort 
has focused on the regimes of prograde orbits with small 
e and I relative to typical WIMPs since most observed 
NEOs have such properties [95-97]. 

However, just like Kozai cycles, secular resonances 
should be able to pull WIMP perihelia outside of the 
Sun if the WIMP orbits originate in the outer layers of 
the Sun, where the orbital precession due to the Sun's 
potential is small. Although there are neither analytic 
nor numerical investigations of secular resonances for 
e > 0.995 relevant for boimd WIMP orbits, extrapolat- 
ing from Williams and Faulkner [98], it appears that for 
fixed a, the prograde resonances lie at higher inclination 
for higher e, so secular resonances will be relevant at high 
inclination, as for Kozai cycles [98]. It is not clear how 
strong these resonances are, although it is unlikely that 
they are much stronger than Kozai resonances. 

Lifetimes: Since Kozai WIMPs dominate the solar- 
captured WIMP DF at the Earth in the simulations, it 
is important to understand the stability of these orbits 
in the true solar system. There are two important ques- 
tions: (i) How long, on average, does it take for a WIMP 
to be perturbed off a Kozai cycle? (ii) How does the 
integrated optical depth per Kozai cycle change? 

Since the diffusion approximation has nothing to say 
about the stability of resonant orbits, we look to simu- 
lations of NEOs again for insight. Unfortunately, NEO 
simulations are either fundamentally short (< 100 Myr) 
or end when NEOs hit the Sun, making it difficult to ex- 
tract estimates of the long-term stability of Kozai cycles. 
There are a few hints from even those short simulations 
with initial conditions significantly different from those 
of WIMPs. First, Gladman et al. [89] find examples of 
NEOs with a < 2 AU in Kozai cycles for tens of Myr in 
their 60 Myr integrations. The lifetimes of those NEOs 
is limited only by the termination of the simulations at 
either 60 Myr or when the body hits the Sun. Thus, it 
seems probable that WIMPs born on Kozai cycles will 
typically stay there for a least of order tens of millions of 
years, and maybe significantly longer. If the timescale to 
perturb a WIMP off a Kozai cycle occurs on timescales 
similar to the ejection timescale (Eq. 104), then WIMPs 
can exist on Kozai cycles of order the age of the solar sys- 
tem. In this case, the DF for WIMPs with a < 1.5 AU 
should be relatively unchanged. On the other hand, if 
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the typical timescale for the removal of a WIMP from a 
Kozai cycle is shorter (such that e/ becomes larger), the 
impact on the DF depends crucially on what timescales 
those WIMPs arc then cither ejected from the solar sys- 
tem or rescattered in the Sun. 

The structural changes to the Kozai cycles in a more 
complex solar system (a is no longer constant, frequent 
switches between librating and circulating modes, Cmax 
and Imax vary) mean that the integrated optical depth 
per Kozai cycles may vary with time (sec Figs. 7 & 8 in 
Gladman et al. [89]). In principle, this could go up or 
down; in the case of the quasi-Kozai cycles in our toy so- 
lar system, the mean optical depth per Kozai cycle went 
up due to occasional periods of very high eccentricity. 
However, given the accessible phase space for WIMPs in 
a more realistic solar system, it is quite possible that the 
mean integrated optical depth per Kozai cycle will go 
down. In this case, the WIMP lifetimes will be length- 
ened, although it is not clear by what amount. 

In summary, we predict that the number density of 
WIMPs with tti < 1.5 AU will be within factors of a few 
of the number densities found in the toy solar system, 
but there are significant error bars in this prediction. We 
find that the additional mechanisms to pull WIMPs out 
of the Sun, angular momentum diffusion or extra secular 
resonances, will at best yield the same number density 
as the WIMPs on Kozai cycles in the toy solar system. 

The total number density will likely depend largely 
on the behavior of WIMPs that were confined to Kozai 
cycles in the toy solar system. The DF will depend 
on the timescales on which WIMPs are removed from 
Kozai cycles, and the timescales for removal from Earth- 
crossing orbits after they have been moved from Kozai 
cycles. If the WIMP-nucleon cross section lies below the 
equilibrium cross section ((7^^ ^ 10^^^ cm^ or a^^ ~ 
10^^" cm^) for the high plateau, perturbations by the in- 
ner planets will reduce the Kozai WIMP DF unless both 
timescales are of order the age of the universe and the 
mean integrated optical depth per Kozai cycle is much 
smaller than in the toy solar system. 

However, for cross sections above the equilibrium cross 
sections, the Kozai WIMP number density will tend to 
increase. If both timescales arc similar to the ejection 
timescale found in Eq. (104), and if the mean opti- 
cal depth per Kozai cycle is similar to what was found 
in the toy solar system, the number density should be 
largely unchanged from what we found in this work. If 
the timescale for rcirnoval of WIMPs from Kozai cycles is 
significantly less than the ejection timescale, or if grav- 
itational perturbations from the planets systematically 
decrease the integrated optical depth per Kozai cycle, 
the DF could be considerably larger. 



2. 1.5 AU< a <2.6 AU 

Given that quasi-Kozai WIMPs have high eccentricity 
and/or high inclination, they also will generically have 



high speed encounters with planets. Thus, we expect the 
timescale for WIMPs to be removed from quasi-Kozai 
orbits, tq, to be similar to that of WIMPs on Kozai cycles. 

After being removed from a quasi-Kozai orbit, the WIMP 
should hit the Sun again in ~ 1 — 10 Myr, according 
to NEO simulations, or get perturbed onto a Jupiter- 
crossing orbit and get ejected. 

We find that e/ ~ tmed/tq, and Et ~ tmax/tmed, where 
tmax is ^0 if the perturbed WIMPs have lifetimes ti > t©, 
and is equal to the median perturbed lifetime otherwise. 
If the other factors in Eq. (84) are ~ 1, this implies that 
any boost or deficit in the number density of WIMPs of 
1.5 AU < a < 2.6 AU goes as n'/n ~ tmax/tq, where n 
is in this case the number density of quasi-Kozai WIMPs 
in the toy solar system simulations. We expect that 
tmax < tQ, and tq > 100 Myr (suggested by the relatively 
short NEO simulations of Gladman et al. [89]). Thus, 
n'/n < 50. If the timescale for the removal of WIMPs 
from quasi-Kozai orbits is greater or the timescale for 
rescattering is less than Iq (quite possible given the fre- 
qTiency with which NEOs hit the Sun), then n'/n will be 
correspondingly less. 



3. Jupiter- Crossing WIMPs 

Before discussing the effects of other planets on 
Jupiter-crossing particles, we summarize the main fea- 
tures of the Jupiter-crossing DF. The plateau in the DF 
was set hy t ^ 10^ yr, with growth in the spikes occur- 
ing at later times due to long-lived Kozai and resonance- 
sticking particles. The vast majority of Jupiter-crossing 
particles are lost by ejection from the solar system rather 
than rescattering in the Sun for a^^ < 10~^^ cm^, al- 
though rescattering becomes more important for larger 
cross sections. 

The outer planets are unlikely to affect the low 
plateau of the Jupiter-crossing DF. Jupiter dominates 
the timescale for Jupiter-crossing WIMPs to be pulled 
out of the Sun; according to Eq. (100), tj, ^ 10^ yr, 
while the timescale for any of the outer planets to re- 
move the perihelion of a passing WIMP is at least an 
order of magnitude longer. Jupiter also has the short- 
est WIMP ejection timescale (Eq. 104) by more than a 
factor of ten. It dominates the Kozai structure of the 
types of orbits on which Jupiter-crossing WIMPs origi- 
nate, and its mean- motion resonances are also by far the 
strongest in the solar system (unless the orbit of the test 
particle is exterior to the orbit of Neptune) [42, 62]. 

However, the outer planets may affect the spikes in 
the Jupiter-crossing WIMP DF because WIMPs that are 
long-lived in the toy solar system may not be long-lived in 
the real solar system. If the orbital node crossings occur 
near one of the outer planets, the WIMP may be quickly 
perturbed from resonant motion and ejected. Thus, it 
is possible that the long-lifetime resonance features in 
the WIMP DF will be less prominent than shown in 
this work, although even in our simulations, the Jupiter- 
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crossing WIMPs are a subdominant contributor to the 
number density. However, shortening the WIMP hfe- 
times in this way only strengthens our conclusion that 
the signal from bound WIMPs in neutrino telescopes is 
unobservably small compared to the signal from unbound 
WIMPs. 

The inner planets may affect the low plateau of the 
bound WIMP DF for the following reason. There is a 
small probability that Jupiter-crossing WIMPs will be 
gravitationally scattered by an inner planet onto an or- 
bit that no longer crosses Jupiter's. If this effect were de- 
scribed by diffusion, the net change to the bound WIMP 
number density would be of order unity; the increase in 
lifetime Et would be canceled by the decrease in et, since 
the timescalc for both scattering in or out of a Jupiter- 
crossing orbit is the same if one inner planet dominates 
the gravitational interaction cross section. However, the 
WIMPs could spread into the low geocentric regions of 
phase space inaccessible in the toy solar system, which is 
important for capture in the Earth. The effect of secu- 
lar or mean-motion resonances on the size of the bound 
WIMP population and the low speed phase space den- 
sity is unclear; resonances could drive WIMPs into the 
Sun, as suggested by NEO and asteroid belt simulations 
[87, 89]. In this case, the importance of resonances de- 
pends on how the typical time for rescattering in the Sun 
relates to the other gravity-dominated times in the solar 
system. 

In siimmary, we siiggest that the height of the low 
plateau of the Jupiter-crossing WIMP DF and the 
Jupiter-crossing WIMP number density will be mostly 
unaffected by the presence of additional planets in the 
solar system, although the inner planets may extend the 
plateau in phase space. We expect that the long-lifetime 
peaks in the Jupiter-crossing WIMP DF will be lower in 
a more realistic solar system due to interactions with the 
outer planets. 



4- Future Simulations 

In order to test the our arguments above and to defini- 
tively determine the bound WIMP DF as a function 
of WIMP parameters, especially at the low geocentric 
speeds inaccessible in the three-body problem but which 
arc so important for WIMP capture in the Earth, we 
would like to perform simulations of WIMP orbits using 
a more realistic model of the solar system. The numerical 
methods presented in Section II and Appendix C should 
be applicable to a more complex solar system with only 
minor tweaking, so we are eager to use our methods for 
future simulations. However, our experience with simu- 
lations in a toy solar system, as well as phenomena high- 
lighted in earlier portions of this section, suggest specific 
challenges to this program. 

The main challenge will be to sample enough orbits to 
have a statistically significant determination of the DF, 
and to do this with finite computational resources. Prom 



our simulations in the toy solar system, wc have learned 
that it is important to determine the long-lifetime tail 
of the WIMP distribution, even if the overall fraction 
of WIMPs in this population is small. The DFs were 
dominated by the small number of particles on Kozai 
cycles (either a < 1.5 AU or on Jupiter-crossing orbits) 
and Jupiter-crossing WIMPs on long-lifetime resonance- 
sticking orbits, about ^ 0.1% of all particles simulated. 
These rare but long-lived WIMPs, especially the Jupiter- 
crossing population, also dominated the uncertainties in 
the DF. However, even getting to this level of uncertainty 
required ^ 10^ CPU-hours per simulation. If we had sim- 
ulated, say, an order of magnitude fewer WIMPs, we may 
not have even identified the long-lived Jupiter-crossing 
population. 

A number of effects we identified earlier in this section 
for orbits in a more realistic solar system will likely affect 
small WIMP populations. For example, the fraction of 
WIMPs with a < 1.5 AU leaving the rescattering peak 
due to angular momentum diffusion will be small: ^ 10~^ 
for a^^ = 10-^^ cm^ and - lO'^ for a^^ = 10'^^ cm^. 
It will be necessary to simulate vast numbers of WIMPs 
with a < 1.5 AU to get good statistics on this population 
and to make sure wc do not miss any important effects. 

We propose the following techniques to maximize the 
statistics on the full solar system bound WIMP DF given 
finite computing time. First, we propose a series of in- 
termediate simulations before simulating WIMPs in com- 
plete solar system to highlight the importance of different 
types of behavior. An initial step may be to simulate or- 
bits in a solar system containing Jupiter, the Earth, and 
Venus (the planets that will likely dominate the behav- 
ior of WIMPs with orbits interior to Jupiter's) on circu- 
lar, coplanar orbits, with the masses of the Earth and 
Venus scaled up by one or two orders of magnitude. We 
choose a low number of planets for simplicity in under- 
standing the simulations, and high masses for the inner 
planets in order to highlight the diffusion processes de- 
scribed in previous sections, for which the gravitational 
cross section scales as Mp (e.g., Eq. 99). The high planet 
masses should shorten the diffusion timescales by a fac- 
tor of Mp'^, which would shorten the total integration 
time. One might then simulate WIMP orbits in a solar 
system with massive inner planets but with more realistic 
planet orbits (highlighting secular resonances), or to sim- 
ulate WIMP orbits in a solar system with the same three 
planets on circular orbits, but for which the masses of 
the Earth and Venus are closer to their true values. One 
could then add the outer planets to the simulation. It 
may be possible to learn how the WIMP DF scales with 
the masses of the inner planets in the simulations with 
higher inner planet masses so the DF could be extrapo- 
lated to small planet masses without needing to simulate 
orbits in a solar system with the true planet masses. Even 
if the latter is not possible, we would learn enough from 
each intermediate simulation to more efficiently run the 
next set of simulations. 

Secondly, we propose weighting the initial conditions 
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to achieve the best statistics with the least amount of 
computational time. The optimal weighting for each in- 
termediate simulation will be guided by the results from 
the previous. For example, say that wc learn from the 
first intermediate stage we propose, a three-planet solar 
system with large inner planet masses, that the popu- 
lation of WIMPs with a < 1.5 AU with perihelia per- 
turbed out of the Sun by angular momentum diffusion 
is significant for neutrino telescope event rates. Perhaps 
the effect is dominated by WIMPs initially scattered into 
a very narrow range of semi-major axes. If we then wish 
to simulate this population in a solar system in which 
the Earth and Venus have their true masses, it makes 
sense to focus the computational resources on this nar- 
row semi-major axis window. Furthermore, in order to 
gain good statistics for this window, we would need at 
least of order 10^ long-lived angular momentum-diffused 
WIMPs. For cr^'^ = 10-''\ e/ - lO"'"^. If we want a sam- 
ple of at least 10'^ WIMPs in this population, we would 
need to simulate ~ 10^ WIMPs. However, ~ 0.1% of 
these WIMPs, or 10'"' total, should initially be on Kozai 
cycles. In order to focus on the angular diffusion popula- 
tion instead of the Kozai population, we would simulate 
all 10^ WIMPs for a short time, 10^ - 10^ years, 
which would be sufficient to identify the Kozai popula- 
tion. At that point, we would only continue simulations 
of the WIMPs not identified as Kozai cycling. Thus, 
we would have good statistics on one WIMP population 
without burning resources on less important populations. 

Therefore, while we believe that getting good statis- 
tics for estimating the event rates in neutrino telescopes 
will be difficult, it will be possible given (i) a clever and 
adaptive simulation strategy, and (ii) patience to acquire 
a sufficient number of CPU cycles. 



C. The Halo Distribution Function 

Throughout the simulations, we assumed that the 
halo WIMP DF was smooth, non-rotating in an iner- 
tial Galactocentric frame (lagging the Sun by a speed 
Vq = 220 km s^^), and had a velocity dispersion of 
a = VqI \pl. These choices are motivated by N-body sim- 
ulations of Milky Way-mass dark matter halos [99, 100]. 
However, there are a few severe limitations to these N- 
body simulations. First, while we hope that the simu- 
lations are a good representation of the real Milky Way, 
there is no way we can directly measure the dark mat- 
ter phase space density. Secondly, these simulations do 
not include baryons, although we know baryons domi- 
nate the gravitational potential within the solar circle. 
Simulations that include a treatment of baryonic disks 
and the accretion of dwarf galaxies suggest that the local 
phase space structure of dark matter depends sensitively 
on the accretion history of the Milky Way [101]. Thirdly, 
dark matter is fundamentally clumpy, with the smallest 
halos corresponding to the size of the free-streaming scale 
[102-104], which for a SUSY WIMP corresponds to about 



M ~ or length scales of ^ 10~^ pc. While high- 
resolution simulations show that very little (~ 0.1 — 0.5% 
[105, 106]) dark matter within the solar circle is in re- 
solved subhalos, these simulations can only probe sub- 
halo masses down to M ~ 10^ - lO^M©. There is thus 
an uncertainty in the degree of dumpiness spanning more 
than 10 orders of magnitude in mass [55]. Here, we de- 
scribe how the DF will change if any of the assumptions 
of our fiducial halo model are challenged. 

We note that the primary change to the DF will be in 
normalization, not shape. The only way to change the 
shape of the bound WIMP DF relative to that calculated 
for our fiducial model for fixed m^^ and elastic scattering 
cross section is to change the distribution of semi-major 
axes or locations of initial scatter in the Sun onto Earth- 
crossing orbits. The former is robust over several orders 
of magnitude in WIMP mass. The latter may be signifi- 
cant for large (m^ > 1 TeV) WIMP masses if the velocity 
phase space is radically different from the fiducial model, 
but will not be significant as long as there is non-trivial 
phase space density of WIMPs at low heliocentric speeds. 

However, the height of the DF is proportional to A^q, 
which is increasingly sensitive to the low speed end of the 
halo DF for increasingly massive WIMPs. This is because 
the halo WIMP energy \% E — (where Voo is the 

heliocentric speed in the absence of the Sun's gravity) but 
the maximum energy a WIMP can lose in a collision with 
a solar nucleus is Qmax = '2n\v^{r)/mA, so it becomes 
hard to scatter high mass WIMPs, high energy WIMPs 
onto bound orbits. If the low speed phase space den- 
sity were increased, Nq would increase, and the bound 
WIMP density would increase relative to the halo den- 
sity. This could be achieved, for example, if the WIMP 
halo were rotating in the same sense as the stellar disk, 
reducing the speed relative speed between the halo and 
the Sun. Conversely, if the low speed halo WIMP density 
were decreased, the bound WIMP population would be 
even more insignificant with respect to the halo. 

While dumpiness in the halo may affect the halo DF at 
the Earth (although it is unlikely that a subhalo is cur- 
rently passing through the solar system [55]), it will have 
surprisingly little effect on the DF of WIMPs bound to 
the solar system if the rate at which clumps pass through 
the solar system is either much higher or lower than the 
equilibrium timescale for the bound WIMP DF. In the 
former case, as long as the velocity distribution of the en- 
semble of subhalos is similar to that of the smooth DM 
component (if the rate at which clumps enter the solar 
system is high) , the bound WIMP DF is proportional to 
the time-averaged capture rate in the Sun, f{v) oc (N^). 
This is unlikely to be significantly different from A^^ cal- 
culated for a purely smooth halo unless the solar system 
is deeply embedded in a dense subhalo. In the latter 
case, passages of a subhalo through the solar system are 
so infrequent that the DF is dominated by the smooth 
component in the halo. 

Diemand et al. [104] estimate that if all Earth-mass 
subhalos survive intact to the present, the rate at which 
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subhalos pass through the solar system is ~ 10~^ yr~^, 
with each passage lasting ^ 50 yr. Diemand et al. [107] 
and Faltenbacher and Diemand [108] find that the veloc- 
ity distribution of subhalos is only slightly biased with 
respect to the smooth component, with the major dis- 
crepancy being a decrement of subhalos with low Galac- 
tocentric speeds due to merging. The escape velocity 
from a subhalo is much smaller than either any charac- 
teristic speed in the solar system or characteristic speeds 
in the solar neighborhood, making it unlikely that the 
Sun is bound to a subhalo. Thus, even if dark matter in 
the solar neighborhood were highly clumpy, the bound 
WIMP DF would resemble that estimated in this work. 



VIII. CONCLUSION 

In conclusion, we highlight the key points of this paper: 

1. We have developed numerical methods to efficiently 
track the highly eccentric solar-captured orbits 
from their initial scatter in the Sun to up to 4.5 Gyr 
without secularly increasing errors in the Jacobi 
constant and without numerical precession. These 
methods will be employed in future simulations of 
WIMPs in a more realistic solar system, and may 
be used to simulate eccentric orbits in other hierar- 
chical systems in which one central body dominates 
the gravitational potential. 

2. We have characterized the bound WIMP DF at the 
Earth as a function of WIMP mass m-^ and spin- 
independent fjp^ and spin-dependent (7^^ elastic 
scattering cross sections. For the range of masses 

= 60 AMU - 500 AMU, we find very lit- 
tle variation in the WIMP DFs aside from the 
mass-dependent rate at which WIMPs scatter onto 
Earth-crossing orbits. In contrast to Damour and 
Krauss [37], we find that the optical depth in the 
Sun to WIMPs imposes a ceiling to the size of 
the WIMP DF. For WIMPs that do not inter- 
sect Jupiter's orbit, the equilibrium DF is reached 
for ~ 10-42 cm2 and ~ lO'^o cm^. 

For WIMPs that intersect Jupiter's orbit, equilib- 
rium is reached for cr^^ ^ 10"^* cm^ or (7^^ ^ 
10-36 cm2. 

3. The maximum phase space density of WIMPs at 
the Earth consistent with current constraints on the 

elastic scattering cross section is significantly less 
than that of WIMPs unbound to the solar system. 
Even though bound WIMPs occupy the low veloc- 
ity phase space that disproportionally contributes 
to the event rates in both direct detection experi- 
ments and neutrino telescopes, the total enhance- 
ment to those event rates is negligible. For direct 
detection experiments, we find that the maximum 
enhancement to dR/dQ occurs at Q = and is 
< 0.5% of the halo event rate. For the XENONIO 



experiment, we predict the maximum enhancement 
integrated over their analysis window is of order 
10~3%. In the MSSM, we find less than order unity 
enhancements to the neutrino-induced muon event 
rate in neutrino telescopes from the annihilation of 
solar-captured WIMPs in the Earth. 

4. Although we only include one planet (Jupiter) in 
our toy solar system, we do not expect that our con- 
clusions would be significantly difi'erent than if we 
had included more planets in our simulations. If the 
other planets are efficient at putting solar-captured 
WIMPs at geocentric speeds v < SO km s~^, there 
may be large increase in the event rate at neutrino 
detectors due to WIMP annihilation in the Earth. 
However, it is unlikely that the boost will be suf- 
ficient to move the event rate above the detection 
threshold for the IceCube neutrino telescope unless 
the halo WIMP DF is significantly different from 
the fiducial model. 

In two other papers in this series, we examine the impact 

of the finite optical depth in the Sun and gravitational 
interactions between WIMPs and Jupiter on the rate of 
WIMP annihilation in the Sun (Paper II); and we char- 
acterize the population of WIMPs bound to the solar 
system by gravitational interactions with Jupiter (Paper 
III). 
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APPENDIX A: WIMP ELASTIC SCATTERING 

1. Spin-Independent Scattering 

For particle physics models of dark matter, the general 
spin-independent ("SI"; scalar) scattering cross section 
has the form [1, 2]: 

^ = ^ [Zf, + iA- Z)Uf FhiQ), (Al) 

where Q is the energy transferred from the WIMP to 
a nucleus of mass vtia (with atomic mass A and charge 
Z) during the scatter, qa is the relative velocity between 
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the particles, fp and /„ arc the proton and neutron ef- 
fective coupHngs to the WIMP, and Fsi{Q) is a nuclear 
form factor. The nuclear form factor used in this set of 
calculations is of the standard exponential form, 



Fsi{Q)=e-Q/'Q\ 
where the coherence energy is 



Qa = 



i.sr 

mAR\' 



(A2) 



(A3) 



and the coherence length (the radius of the nucleus A) is 
set to 

Ra = 1 fm[0.3 + 0.91(mA/(GeV/c2))V3]. (A4) 

The nuclear form factor quantifies the extent to which the 
WIMP interacts coherently with the nucleus as a whole 
(if the de Broglie wavelength of the nucleus is small), or 
incoherently with the nucleons individually. 

It is often more convenient to use the center-of-mass 
differential cross section. Using the functional form of 
the energy transfer 



rriA 



cos 6 



(A5) 



where 
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the differential cross section is 

da^^ _ 1 dQ da 
dn ~ 2^ d{cos 6) dQ 
_ 1 iJi\ 2 ( do- 
2^^^^ [dQ 

47r TT 

^ ■F^(Q(cosg)) 

47r 



(A6) 

(A7) 
(A8) 

Z)UfF\Q) (A9) 
(AlO) 



We have parameterized the strength of the interaction by 
OA- If fp = fn, which is often a good approximation for 
both supersymmetric and UED models. 



^SI " ,;a \'i fi 

"A — ~I^A^ Jm 
TT 



(All) 



so that the strength of the coupling between a nucleus 
and the WIMP depends only on the atomic number of 

the nucleus. This coupling can also be parameterized in 
terms of the strength of the WIMP-proton (or -neutron) 
cross section: 



SI _ Ma a2^SI 



(A12) 



which is useful since experimental constraints on the spin- 
independent cross section are reported in terms of the 



WIMP-nucleon cross section. In the limit of high WIMP 
mass, 



Ma 
Mp 

<^A 



ruA 



ml P 
A'a'p', 



(A13) 
(A14) 

(A15) 
(A16) 



where the last approximation can be made since tua 
Am„. 



2. Spin-Dependent Scattering 

The likely WIMP candidates for both the MSSM (neu- 
trino x) a-iid UED (Kaluza-Klein photon B^^'') theories 
can have elastic axial- vector interactions with quarks, via 
squarks in the MSSM or the lightest Kaluza-Klein exci- 
tation of quarks q^^^ in UED models. In both cases, the 
spin-dependent (SD) WIMP interaction with a nucleus 
of atomic number A can be parameterized as [1, 3] 



dQ 



= a X 



2mA 



AV(J+l)f|^(|q|: 



(A17) 



where 




TT? — 2 1 — \^ 



MSSM 
UED 



(A18) 



parameterizes the coupling in each theory. Here, g' is the 
coupling constant for the B boson in electroweak theory, 
and m^ci) and m^(i) are the masses of the B^^^ and g^^^ 
particles respectively. The other quantities in Eq. (A17) 
depend on nuclear properties. Here J is the total angular 
momentum of the nucleus, and 



A = J [ap(S'p) -h a„(S'„)] , 



(A19) 



where a„ and describe the WIMP couplings to the neu- 
tron and proton, and (S^) and (S'p) are the spin expecta- 
tion values for the neutrons and protons within the nu- 
cleus. The couplings a„ and are derived from specific 
WIMP models, while the spin expectation values must 
be calculated using detailed nuclear physics models [e.g., 
1, 109-111], and calculations using different techniques 
often yield different results. The function F5£)(|q|) is the 
spin-dependent nuclear form factor as a function of the 
momentum transfer |q|. Its form must be carefully cal- 
culated for each nucleus of interest [112, and references 
therein] . 

There are several important differences between the 
form of the spin-dependent and spin-independent cross 
sections that have major implications for detection ex- 
periment design. The first point is that nuclei with even 



38 



numbers of protons and neutrons will have zero spin- 
dependent interactions with WIMPs. Secondly, the spin- 
dependent cross section has a much weaker dependence 
on the atomic mass than the spin-independent cross sec- 
tion. This is apparent if Eq. (A17) is written in the same 
form as Eq. (AlO), 

da^^ 1 dQ da^^ 



dO 



27rdcos6» dQ 

2iT m,A 



„ J( J + l)aA^ 



4tt 



where 



= -ij\j{J+l)aA\ 
In the limit that mwiMP ^ '^A, 

„SD ^ a2 



unlike 



(A20) 

(A21) 
(A22) 

(A23) 

(A24) 

(A25) 

for the spin-independent case. Therefore, even if cTp ^ > 
CTp-^ or a^^ > crf^, the spin-independent cross section 
may dominate for heavy nuclei. The spin-dependent 
cross section could be large if J scaled with A (since 
a A oc J^), but this is not the case for heavy nuclei. Note 
that, in contrast to predictions for spin-independent scat- 
tering, the spin-dependent WIMP-proton and WIMP- 
ncutron cross sections are generally not the same to 
within a few percent. 



APPENDIX B: SUBSEQUENT SCATTERING IN 
THE SUN 



Each time a particle passes through the Sun, there is 
a probability 



af oc A" 



1 



(Bl) 



that it will be scattered at least once, given the optical 
depth r for one jaunt through the Sun. Since the WIMP- 

nuclcon cross sections relevant to this paper imply low 
opacity in the Sun (r ^ 10~^), the scattering probability 
per solar passage is well approximated by 



Pscatt — 1 ~ (l 

Ri r. 



(B2) 
(B3) 



Instead of calculating the scattering probability r on 
the fly, we create a table for optical depth indexed by the 
semi-major axis and Kepler perihelion of the orbit, and 



then interpolate for a particular orbit through the Sun. 
The optical depth in differential form is given by 



dr 
d/dQ 



E 



dTA 

dldQ 



daA 
dQ ' 



(B4) 
(B5) 



where / denotes the particle trajectory, UAil) is the num- 
ber density of element A in the Sun at position I along 
the path, and dcr^/dQ is the differential elastic scattering 
cross section with respect to the energy transfer Q be- 
tween element A and the WIMP. Since we assume that 
spin-independent scattering dominates in the Sun, the 
integral over energy transfer can be computed using the 
form of the differential cross section in Eq. (Al) and the 
form factor in Eq. (A2): 

^ (Be) 



dr 
d7 



A 



'TTv{iy 



{A - Z)UfQA (B7) 



i/QA 



where wc have used the approximation of a zero-temper- 
ature Sun to set Vrei = v{l). Using Eq. (62), we find the 
maximum energy transfer 



Omax,A = 2^v{lf 

ruA 



(B8) 



The integration of Eq. (B8) is greatly simplified be- 
cause the torque on the particle by Jupiter is negligible 
in the Sun compared to the rest of the orbit. Therefore, 



dl = v{t)dt 

= vit{r)) 
v{r{t)) 



dt 
dr 



dr 



\vr{rm 



dr, 



where 



v{E,r) = ^/2[E-^Q{r)] 
is the particle's speed and 

\vriE,J,r)\ = ^2[E-^Qir)]-jyr^ 
is the radial velocity of the particle. Thus, 

dT{E,J) _ v{r) dr 
dr br(r)| dl 

and the total optical depth along the path is 
4 

TT 



(B9) 
(BIO) 

(Bll) 

(B12) 
(B13) 

(B14) 



TiE, J) = - V mAQA[Zfp + {A- Z)Uf 
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Re riAir) ( 1 - e-2Mi'''(s,r-)/m^QA\ 
dr- ^ ^ 
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(B15) 
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In order to express the optical depth t as a function of 
tie semi-major axis and eccentricity, We use the relations 



E = ±^ 
2a 

J2 = ±GMea(e2 



1), 



(B16) 
(B17) 



where Mc — hMq is the central mass, as determined by 
Eq. (34), and the upper (lower) sign is used for hyper- 
bolic (elliptical) orbits. Therefore, 



v{a,r) = v{a,r)/\/GMc 



(B18) 



\v{a,e,r)\ = \vr{a,e,r)\/^/GM^ 



where $0 = ^q/GMc- If we insert these expression into 
Eq. (B15), 

r(a,e) = iJ—J2^AQA[Zfp + iA-Z)U]^ (B20) 



X / dr — 



(r) (1 



-2iiAGM^v^{a,r)/mAQA 



) 



v{a,e,r)\vr{a,e,r) 



We make a look-up table for r using for the choice /i = 1 . 
and then scale r by a factor of /U"^. There is also a 
factor of n in the exponent. However, its impact on r is 
negligible since |/i — 1| ^ 10^^ — 10^'''. 

If the particle scatters in the Sun, its new phase space 
coordinates can be determined by sampling the scattering 
distribution 



drdn 



\vr{E,J,r)\ dn 



where O is the center-of-mass scattering solid angle. 



APPENDIX C: DISTRIBUTION FUNCTION 
ESTIMATORS 

In this section, we describe the outputs of the simula- 
tions, and how to estimate the bound distribution func- 
tion from these data. 

Our method is to find the average DF along Earth's 
path. We record the phase space coordinates of particles 
passing near the Earth's orbit. Since we treat the Earth's 
orbit as circular and coplanar with Jupiter's orbit, this 
means that we focus on particles passing through the wall 
of a cylinder of height 2zc centered on the reference plane 
and radius from the Sun. Thus, the raw data product 



is the flux of dark matter particles through the Earth's 
orbit as a function of time. 

To convert the flux at position x and time i, F{x,t), 
into a DF /(x, v, t), we assume that the timcscalc of vari- 
ation in the distribution function is much larger than the 
typical dynamical timescale of particles in the solar sys- 
tem (^ year). We adopt the usual argument [cf. 113] 
to relate the flux as a function of velocity dF/dv to the 
distribution function. Consider particles passing outward 
through a wall of area 5 A with a unit vector normal to 
the surface n. For particles with velocity between v and 
V -I- (5v, the particles that pass through the wall in time 
St inhabit a prism volume of base 6 A, long side v5t, and 
height Stv-n. The total number of particles with velocity 
between v and v + (5v passing out through the surface 
SA per unit time 6t is 



dF(x, t) 
dv 



dvSA6t = f{x,v,t){v6t) ■ {5Ah)dv (CI) 
= /(x, V, t)v cos -fdvSA6t, (C2) 



where cos 7 = v • ii/v. In the simulations, we do not care 

if the particles pass inward or outward through the wall 
of the cylinder, so we estimate the distribution function 
from the simulations using 



di^(x, t) 



dv 



dvSA6t = / (x, V, t) V I cos 7 1 dv6A6t, (C3) 



or 



/(x,v,i) 



dF(x,t)/dv 



V cos 7 
dF(x, t)/dv 



(C4) 
(C5) 



since Vr = t;cos7 is the velocity component normal to 
the wall of the cylinder (i.e., the radial component of the 
velocity) . 

We now describe in detail how to estimate the distribu- 
tion function from the data obtained in the simulations. 
For each simulation, we start integrating the orbits of Np 
particles (Table I) at time ti since the birth of the solar 
system. Particles scatter onto bound. Earth-crossing or- 
bits at a rate NQ{ti), where ti is the time at which the 
particle first scatters onto a bound orbit. In principle, 
Nq can vary with time if the halo dark matter distribu- 
tion function varies on timescales shorter than the age of 
the solar system, but we assume that the halo distribu- 
tion function is static, so that NQ{ti) = Nq^. 

Each time a particle a crosses through the cylinder 
wall, we record the time of passage ta^ (here, /3 labels the 
particular passage of the particle a through the Earth's 
orbit) since the start of the simulation at ti, position x^,^, 
and velocity Vap- The height Zc is chosen to be larger 
than the radius of the Earth in order to improve 
statistics, but is small enough {zc «C 1 AU) so that the 
estimate should be unaffected by gradients in flux as a 
function of height above the reference plane. 
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Each particle crossing can be characterized as one 
point in a six- dimensional phase space: ria/j, the vec- 
tor describing the orientation (0, z) of the particle when 

it crosses the cylinder of radius a® ; the three components 
of the velocity v^/j; and ta0- The vector Hq,^ only has 
two independent coordinates since the radial component 
of Xq,/3 is fixed. We estimate the flux of particles passing 
passing through a patch of the cylinder at position n in 
the cylinder at time t since the birth of the solar system, 
for which the particles had initial scattering time in the 
Sun at time U, with velocity between v and v + dv, as 



the flux is measured, 



dF 
dvdti 



1/ / dX^w{X)S{X- Xa) 



a=l/3=l 

y -•Va0,t- (ti+tafs)) (C6) 

for each experiment. Here, F denotes that this is an 
estimator for the true flux F. The total flux can be es- 
timated by integrating Eq. (C6) over ti and v. Na is 
the total number of times particle a crosses the Earth's 
orbit. The weight function w{X) describes how we sam- 
ple the initial conditions A relative to the initial particle 
distribution at the first scatter. The denominator of Eq. 
(C6) normahzes the flux. 

Since we sample the bound. Earth-crossing WIMPs to 
the same density as they scatter onto such orbits in the 
solar system, w = 1 for each particle a. Thus, 



^w{X)S{X-Xc,) = ^^(A-A„) (C7) 



a=l 



SO that 



dXj2 wiX)SiX - A«) = Np, 



(C8) 



where the integral over A spans the entire range of Aa- 
The flux at position n as a function of velocity, observa- 
tion time, and initial time ti is 



dF N '^'^ 

X] X] ^^'^^ (n - n„/3, V - 



^ a=l /3=1 



t-{ti+t„0)). (C9) 



We are interested in the flux arising from particles en- 
tering the solar system at all times prior to the present, 
not just at a particular time ti. Therefore, to estimate 
the total flux in a unit volume of velocity-space, one must 
integrate Eq. (C6) over ti, in the range between the time 
of the formation of the solar system and the time at which 



dF 
d^ 



-L 



, dF 
dvdij 



(CIO) 



^ ^ (5^^^ (n - n„/3, V - v«/3) 



Nr, 

X Q{t-tcl3) 

In order to get better statistics for the flux through the 
Earth, we average the flux in Eq. (CIO) over all positions 
n on the cylinder wall. In this case. 



d n = = 2 X 2-Kai^Zc, 



(Cll) 



cylinder 



the whole area through which we count particle crossings. 
This implies that the averaged flux is 



dF(n, t) _ J_ 

5A 



dv 



d^n 



cylinder 



dF 
d7 



(C12) 



N 1 ^° 

P a = l/3=l 

X Q{t-t„0). 



(C13) 



In effect, we are averaging the flux over the Earth's orbit. 
We flnd the local estimate of the distribution function by 
inserting Eq. (CI 3) into Eq. (C5). 

To find the distribution function in the frame of the 
Earth, we make a Galilean transformation u = v — v®, 
where v^ is the circular velocity of the Earth about the 
Sun, to flnd 



x,u,t) = /(x,u-|- ve,f). 



(C14) 



1. Estimating Distribution Functions in Practice 

In practice, there are 10^ — 10^ Earth-orbit crossings 
in each simulation. In order to present and use the DFs 
in a manageable form, we use a small Zc and bin the dis- 
tribution function in velocity space. We set Zc = lOiJ®, 
but using different Zc up to Zc = 10~^ AU (the largest 
value we tried) yields consistent DFs, demonstrating the 
desired result that the estimate for the DF does not de- 
pend on the choice of Zc- 

The most straightforward way of estimating uncer- 
tainty in the distribution function and any calculations 
derived from it is to use bootstrap resampling. Bootstrap 
resampling yields accurate parameter and error estima- 
tion if the data sample the underlying distribution well. 
In each resampling, we select Np initial conditions with 
replacement from the Np WIMPs. We then calculate all 
distribution functions and event rates using the trajecto- 
ries and crossings of the new sample as described in the 
previous section. 
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2. The Distribution Function in the Earth 

In the previous section, we found DFs in the absence of 
the Earth's gravity. However, since both direct detection 
experiments and neutrino telescopes are sensitive to par- 
ticles within the potential well of the Earth, it is neces- 
sary to find the mapping between the velocity coordinates 
at distances <C 1 AU from the Earth but well outside the 
Earth's gravitational field and those at distances at which 
the Earth's gravity is significant. Let v = (w, 0, (f>) denote 
the velocity outside the Earth's gravitational field in an 
inertial frame centered on and moving with the Earth, 
with the polar axis along the Earth's direction of mo- 
tion, and the velocity v/oc = {vioc,0ioc,4>ioc) be in the 
Earth's gravitational field at a position R = {R, ip) 
from the Earth's center, where v/oc is also in an inertial 
frame centered on and moving with the Earth. In these 
coordinates, the angles 9, 9ioc, and C are measured rela- 
tive to the direction of motion of the Earth with respect 
to the Sun, and the (f), (f)ioc, and tp angles are azimuthal. 

Since the particle energy E and angular momentum J 
with respect to the Earth are approximately conserved 
near the Earth, the local DF /loc of dark matter in the 
gravitational field of the Earth can be written as 

/;oc(R, ^loc) = /(v(Vioe, R)). (C15) 

Here, /(v) is the dark matter DF in the frame of the 

Earth but far from the Earth's center. Eq. (C15) is a 
restatement of Liouville's theorem. The number of parti- 
cles in an interval between (R, vioc) and (R -|- dR, v;oc + 
d-vioc) is 

dN = fioc(R, Vloc)<i'lid\loc. (C16) 

If the DF /(v) were isotropic, then the mapping be- 
tween velocity coordinates would be greatly simplified. 
In such a situation, the speeds v and vioc are related 
through conservation of energy, 

E=lv' = lvUR)+^M, (C17) 

assuming that the Earth's potential $® is spherical. 
Therefore, the number of dark matter particles with po- 



sitions between R and R + dR and speeds between v/oc 
and vioc + dvioc would be represented as 

dNi,o = 4TrvlJiv{R, vioc))d'lidvioc. (C18) 

However, the DFs are not isotropic in the frame of the 
Earth. Thus, it is necessary to find v in terms of the 
velocity v;oc at position R. The speeds are still related 
by Eq. (C17) , so that i; is a function of only two variables, 
vioc and R. The angular coordinates {9, (j)), however, will 
now be a complicated function of all six local phase space 
coordinates, so that the number of particles at (R, v/oc) 
is described as 



dN = f{v{R, vioc),9{R, v,o,), 0(R, Vioc)) 

X R'^vfoc<iR<icos(diljdviocd cos 9iocd(pioc- (C19) 

To relate the angular coordinates, we make use of an- 
gular momentum conservation as well as energy conser- 
vation, and the fact that the problem reduces to a spheri- 
cally symmetric two-body problem. Since orbits are con- 
fined to a plane, R and vioc are a set of basis vectors for 
the orbital plane if the vectors are not parallel. Then, 
in general, the position R/ar and velocity v far from the 
Earth can be described by 

Rfar = aR + f3vioc, (C20) 
V = jR + Svioc, (C21) 

where the coefficients a, (3, 7, and 5 only depend on the 
local coordinates R and vjoc, and J. If the Earth's 
potential were purely Keplerian, a and (3 would be the 
Gauss / and g functions [see Section 2.5 in 57], with 
7 = d and S = $. The functional form of the coefii- 
cients is different in the case of non-Keplerian spherically 
symmetric potentials, but the general framework of Eqs. 
(C20) and (C21) holds. Therefore, Eqs. (C20) and (C21) 
describe the mapping between coordinates in the gravi- 
tational field of the Earth to those outside the Earth's 
sphere of influence. 
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